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Abstract 


This  paper  presents  numerical  results  for  a  large  and  varied  set  of  problems  using 
software  that  is  widely  available  and  has  undergone  extensive  testing.  The  algorithms 
implemented  in  this  software  include  Newton-based  linesearch  and  trust-region  methods 
for  unconstrained  optimization,  as  well  as  Gauss-Newton,  Levenberg-Marquardt,  and 
special  quasi-Newton  methods  for  nonlinear  least  squares.  Rather  than  give  a  critical 
assessment  of  the  software  itself,  our  original  purpose  was  to  use  the  best  available 
software  to  compare  the  underlying  algorithms,  to  identify  classes  of  problems  for  each 
method  on  which  the  performance  is  either  very  good  or  very  poor,  and  to  provide 
benchmarks  for  future  work  in  nonlinear  least  squares  and  unconstrained  optimization. 
The  variability  in  the  results  made  it  impossible  to  meet  either  of  the  first  two  goals; 
however  the  results  are  significant  as  a  step  toward  explaining  why  these  aims  are  so 
difficult  to  achieve.  (  (Cp-  )  / _ 
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1.  Introduction 


% 


The  nonlinear  least-squares  problem  is  that  of  minimizing  a  sum  of  squares 

s".  \ 

i=l 

in  which  each  fa  is  a  real-valued  function  having  continuous  second  partial  derivatives. 


The  problem  can  also  be  posed  as  a  minimization  of  the  I2  norm  of  a  multivariate 


function : 


min  \  ||/(*)||5, 


where 


/(*)  = 


( <h(x)  \ 

<h(x) 


We  shall  refer  to  the  function  ^  ||/(z)||2  as  the  nonlinear  least-squares  objective  func¬ 
tion.  It  is  assumed  that  n  and  m  are  relatively  small,  so  that  algorithms  are  not  for¬ 
mulated  with  any  special  considerations  for  limited  storage.  A  common  instance  is  the 
choice  of  parameters  within  a  nonlinear  model  <p : 


«=1 

where  d,  are  observations  at  prescribed  values  T{. 


Many  specialized  algorithms  have  been  developed  to  take  advantage  of  the  structure 
of  the  nonlinear  least-squares  objective.  This  paper  presents  numerical  results  for  a  large 
and  varied  set  of  problems  using  "'ftware  that  is  widely  available  and  has  undergone 
extensive  testing.  The  algorithms  implemented  in  this  software  include  Newton-based 
linesearch  and  trust-region  methods  for  unconstrained  optimization,  as  well  as  Gauss- 
Newton,  Levenberg-Marquardt,  and  special  quasi-Newton  methods  for  nonlinear  least 
squares.  Rather  than  give  a  critical  assessment  of  the  software  itself,  our  original  purpose 
was  to  use  the  best  available  software  to  compare  the  underlying  algorithms,  to  identify 
classes  of  problems  on  which  the  performance  of  each  method  is  either  very  good  or 
very  poor,  and  to  provide  benchmarks  for  future  work  in  nonlinear  least  squares  and 
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unconstrained  optimization.  The  variability  in  the  results  made  it  impossible  to  meet 
either  of  the  first  two  goals.  However,  the  results  are  significant  in  that  they  reveal 
a  great  deal  about  the  reasons  these  aims  why  aims  are  so  difficult  to  accomplish. 
For  treatment  of  issues  and  methodology  for  software  performance  evaluation  see,  e.  g., 
More,  Garbow,  and  Hillstrom  [1978;  1981],  Hiebert  [1979;  1981],  and  Hanson  and  Krogh 
[1987].  Hiebert  [1979;  1981]  conducts  an  extensive  evaluation  of  twelve  programs  for 
nonlinear  least  squares,  in  which  she  includes  software  that  uses  first  derivatives  as  well 
as  some  that  does  not.  In  the  present  study,  all  of  the  software  requires  first  if  not 
second  derivatives  of  the  problem  functions. 

This  paper  is  organized  as  follows.  Section  2  reviews  computational  techniques 
for  the  unconstrained  optimization  problem.  These  methods  are  of  interest  because 
the  nonlinear  least  squares  problem  is  a  particular  instance  of  unconstrained  optimiza¬ 
tion,  so  that  special-purpose  algorithms  for  sums  of  squares  should  compare  favorably 
in  performance  with  those  developed  for  the  more  general  case.  Moreover,  much  of  the 
motivation  for  unconstrained  optimization  methods  is  also  relevant  to  algorithm  devel¬ 
opment  for  nonlinear  least  squares.  Our  emphasis  is  on  computational  issues  related 
to  the  methods  included  in  this  study.  For  more  extensive  treatment  of  unconstrained 
optimization  algorithms,  see  Fletcher  [1980],  Gill,  Murray,  and  Wright  [1981],  Dennis 
and  Schnabel  [1983],  and  More  and  Sorensen  [1984],  Section  3  reviews  the  principal 
approaches  that  are  used  in  software  for  nonlinear  least-squares  problems.  These  are 
Gauss-Newton  methods;  Levenberg-Marquardt  methods,  one  of  which  is  implemented  in 
the  software  package  MINPACK  [More  (1978),  More,  Garbow,  and  Hillstrom  (1980)]; 
corrected  Gauss-Newton  methods  [Gill  and  Murray  (1978)],  which  are  the  basis  for  the 
NAG  Library  nonlinear  least-squares  software;  and  methods  that  form  quasi-Newton 
approximations  to  the  term  B  =  in  the  nonlinear  least-squares  Hessian, 

a  strategy  that  is  adaptively  combined  with  a  Gauss-Newton  method  and  a  Levenberg- 
Marquardt  method  in  the  computer  algorithm  NL2S0L  [Dennis,  Gay,  and  Welsch  (1981a, 
1981b)].  We  assume  a  knowledge  of  numerical  techniques  for  linear  least-squares  (e.  g., 
Lawson  and  Hanson  [1974],  and  Golub  and  Van  Loan  [1983]).  For  more  information 
on  algorithms  specific  to  nonlinear  least-squares  problems,  see  Fraley  [1988]  and  the 
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references  cited  in  that  paper.  Section  4  is  a  summary  and  discussion  of  the  numerical 
results.  Section  5  contains  tables  of  all  of  the  results,  as  weH  as  information  about  the 
software  and  test  problems  used  in  obtaining  them. 
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1.1  Definitions  and  Notation 

We  shall  use  the  following  definitions  and  notational  conventions  : 


•  Generally  subscripts  on  a  function  mean  that  the  function  is  evaluated  at  the 
corresponding  subscripted  variable  (for  example,  /*  =  /(x*)).  An  exception  is 
made  for  the  residual  functions  ij>i,  where  the  subscript  is  the  component  index  for 
the  vector  /. 

•  /  -  The  vector  of  nonlinear  functions  whose  I2  norm  is  to  be  minimized. 

The  nonlinear  least-squares  problem  is 

jsi?.  5 

where  the  factor  |  is  introduced  in  order  to  avoid  a  factor  of  two  in  the  derivatives. 


•  fa  -  The  ith  residual  function,  also  the  ith  component  of  the  vector  /. 

f  Mx) 

/(*)  =  ; 

An  alternative  formulation  of  the  nonlinear  least-squares  problem  is 

5  X>M!. 

*=1 

where  each  <f>i(x)  is  a  smooth  function  mapping  SJ”  to  JJ. 


•  J  -  The  Jacobian  matrix  of  /. 

J(x)  =  Vf(x)  = 


9<t>  1 

9<t>  1 

9ij 

dx„ 

d*  m 

Mm. 

dx\ 

9xn 

•  g  -  The  gradient  of  the  nonlinear  least-squares  objective. 

g(x)  =  V  Q/(x)T/(z)^  =  J(x)rf(x) 
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•  B  -  The  part  of  the  Hessian  matrix  of  the  nonlinear  least-squares  objective  that 
involves  second  derivatives  of  the  residual  functions. 

Q/MT/(*))  =  +  B(x), 

where 

m 

B(x)  = 

t=i 

•  ^(A)  -  The  range  of  A. 

If  A  is  an  m  x  n  matrix,  then 

H(A)  =  {b  \  Ax  =  b  for  some  x  G  5?”} 
is  a  subspace  of  9?m. 

•  M{A)  -  The  null  space  of  A. 

If  A  is  an  m  x  n  matrix,  then 

Af(A)  =  {z  |  Az  =  0} 

is  a  subspace  of  Rn.  A'X.d)  is  the  orthogonal  complement  of  7l(AT)  in  S?n. 
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2.  Methods  for  Unconstrained  Optimization 

2.1  Local  Quadratic  Approximation 

Software  for  the  unconstrained  optimization  problem 

min  /Yx), 

is  usually  based  on  successive  minimization  of  a  quadratic  approximation 

Qk(p)  =  V?Jp+\pTHkp  (2.1.1) 

for  T{xk  +  p)  -  ?{xk),  the  change  in  T  at  xk.  The  matrix  Hk  is  positive  definite,  so 
that  the  model  Qk(p)  has  a  well-defined  minimum 

Pfc  =  -n^vrk, 


that  can  be  computed  efficiently.  Positive  definiteness  of  Hk  also  means  that  pk  is  a 
descent  direction  for  T  at  xk,  which  is  essential  for  linesearch  methods  (see  Section 
2.2.1).  For  methods  based  on  (2.1.1),  the  condition 


lim 

Ar— »oo 


|(gfc-V^(x»))Pfc| 

l|P*ll 


=  0 


(2.1.2) 


is  equivalent  to  local  superlinear  convergence  of  the  sequence  {xk  +  pk}  to  an  isolated 
local  minimum  x *  of  T  (see  Dennis  and  More  [1974;  1977]).  The  relationship  (2.1.2) 
implies  that  the  step  pk  approaches  the  Newton  step  in  both  magnitude  and  direction, 
although  the  sequence  of  matrices  {Hk}  need  not  converge  to  V2T{x*). 


Section  2.2  is  concerned  with  modifications  that  are  used  to  enforce  convergence 
from  an  arbitrary  starting  point.  These  modifications  fall  into  two  categories:  linesearch 
methods  and  trust-region  methods.  Section  2.3  deals  with  the  choice  of  Hk  in  (2.1.1), 
so  that  condition  (2.1.2)  for  superlinear  convergence  is  satisfied.  We  discuss  methods 
that  use  exact  second  derivatives  as  well  as  quasi- Newton  approximations. 


2.2  Globalization  Strategies 
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2.2.1  Linesearch  Approach 


Linesearch  methods  obtain  a  new  iterate  in  two  essentially  separate  phases.  First, 
a  descent  direction  pk  is  found  for  T\  that  is,  a  vector  pk  is  computed  for  which 

VJFjp*  <  0.  (2.2.1) 

Condition  (2.2.1)  is  equivalent  to  saying  that  T  initially  decreases  along  the  direction 
Pk  from  x k-  Various  ways  of  defining  pk  are  discussed  in  Section  2.3.  This  section  is 
concerned  with  the  second  phase  of  a  linesearch  method,  that  of  finding  a  steplength 
ajt  satisfying 

T(xk  +  QkPk)  <  ?{xk),  (2.2.2) 

once  a  descent  direction  is  obtained. 

Because  of  (2.2.1),  condition  (2.2.2)  can  be  satisfied  by  choosing  a  sufficiently 
small  value  of  a*,  but  the  result  may  not  be  an  appreciable  reduction  in  T ’.  In  fact, 
{a;*  +  OfcPfc}  may  converge  to  a  point  that  is  not  a  stationary  point  unless  conditions 
stronger  than  (2.2.2)  are  imposed  on  a*  (see,  e.  g.,  Dennis  and  Schnabel  [1983],  Chapter 
6).  On  the  other  hand,  finding  a  minimum  of  T  along  pk  is  an  iterative  process  which 
could  require  many  function  and  derivative  evaluations.  Steplength  algorithms  instead 
compute  Qk  that  satisfies  conditions  sufficient  to  ensure  convergence  to  a  stationary 
point  of  T  whenever  the  sequence  {pk}  is  bounded  away  from  orthogonality  to  the 
gradient. 

The  work  of  Goldstein  [1965;  1967],  Armijo  [1966],  Goldstein  and  Price  [1967],  and 
Wolfe  [1969;  1971],  (see  also  Ortega  and  Rheinboldt  [1970])  established  the  fundamen¬ 
tal  principles  underlying  most  steplength  algorithms.  A  simple  strategy  for  sufficient 
decrease  is  based  on  the  condition 

F(xk  +  <*kPk )  ~  H*k)  <  pak'VT'kPk,  (2.2.3) 

for  p  €  [0,1).  An  initial  value  (usually  a*  =  1)  is  tried  first,  and  then  a  backtracking 
strategy  is  used  to  reduce  it  until  an  admissible  step  is  found.  The  steplength  strategy 
of  Gill  et  al.  [1979],  combines  (2.2.3)  with  the  condition 

\VT(xk  +  OfcP*)TPk|  <  -T)VfJpk,  (2.2.4) 
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for  77  €  [0,1),  which  keeps  the  steplength  bounded  away  from  zero  by  forcing  it  to 
approximate  a  local  minimum  of  T  along  pk ■  A  procedure  for  one-dimensional  mini¬ 
mization  is  truncated,  using  (2.2.4)  as  the  criterion  for  termination.  This  is  accomplished 
by  polynomial  interpolation  to  the  function 

$(a)  =  T{xk  +  apk),  (2.2.5) 


together  with  some  safeguards  to  prevent  iterates  from  being  either  too  close  together 
or  too  far  apart.  An  exact  minimization  would  be  carried  out  for  tj  =  0  in  (2.2.4), 
while  larger  values  of  77  increasingly  relax  this  requirement.  When  p  <  77,  an  interval 
of  steplengths  satisfies  both  (2.2.3)  and  (2.2.4);  if  p  is  chosen  sufficiently  small,  then 
(2.2.3)  almost  always  holds  when  (2.2.4)  does.  When  p  >  tj,  a  backtracking  strategy 
may  be  used  if  (2.2.3)  fails  to  hold  for  the  steplength  computed  in  the  one-dimensional 
minimization.  If  <  0  and  a*  satisfies  (2.2.3)  and  (2.2.4),  then 


lim 

fc— ►  00 


V^TPfc 

W2 


=  0, 


which  implies  convergence  to  a  stationary  point  of  JF  provided  { pk }  remains  uniformly 
bounded  away  from  orthogonality  to  {V/*}.  If  p  <  0.5,  both  conditions  (2.2.3)  and 
(2.2.4)  are  automatically  satisfied  by  superlinearly  or  quadratically  convergent  algorithms 
with  a*  =  1  when  xk  is  sufficiently  close  to  a  local  minimum. 

Although  the  theory  allows  considerable  flexibility  in  choosing  the  interpolant  to 
$(a)  and  other  parameters  in  the  univariate  minimization,  as  well  as  in  the  choice  of  p 
and  77  in  (2.2.3)  and  (2.2.4),  in  practice  performance  on  difficult  problems  may  be  very 
sensitive  to  these  factors.  Moreover,  safeguarding  in  univariate  minimization  requires 
specification  of  a  finite  interval  of  uncertainty  in  which  the  minimum  is  presumed  to  lie. 
If  pk  is  very  large,  it  could  happen  that  no  satisfactory  approximation  to  a  minimum 
along  that  direction  can  be  found,  resulting  in  an  excessively  small  step. 


2.2.2  Trust-Region  Approach 

Trust-region  methods  were  first  developed  for  nonlinear  least  squares  [Levenberg 
(1944);  Morrison  (1960);  Marquardt  (1963))  (see  Section  3.2),  and  later  independently 
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for  general  unconstrained  minimization  [Goldfeld,  Quandt,  and  Trotter  (1966)].  Moti¬ 
vation  for  trust-region  methods  comes  from  the  following  observation:  if  the  step  to  the 
unconstrained  minimum  of  the  current  local  model  for  T(x  +  p)  —  T(x)  is  relatively 
large,  then  it  probably  falls  outside  the  region  in  which  the  model  is  applicable.  The 
basic  idea  is  to  define  a  neighborhood  of  the  current  point  over  which  an  approximate 
minimization  of  a  local  model  of  the  change  in  T  will  yield  a  suitable  step  to  the  next 
iterate. 

The  local  model  and  constraints  defining  the  neighborhood  are  chosen  so  that  the 
subproblem  has  a  well-defined  minimum,  and  so  that  fast  local  convergence  is  possi¬ 
ble  with  the  unconstrained  model.  Typically,  the  model  at  Xk  is  a  quadratic  function 
V-Fjp-F  |  pT Hkp,  and  an  upper  bound  is  imposed  on  a  scaled  I2  norm  of  p,  giving  the 
subproblem 

mmVfJp+±pTffkp  (2.2.6) 

subject  to||Z)jfcp||2  <  6k. 

For  practical  reasons,  the  scaling  matrices  Dk  are  usually  diagonal  (with  positive  diagonal 
entries).  Solving  (2.2.6)  is  equivalent  to  minimizing  the  quadratic  function 

V^Jp+  \pT  ( Hk  +  X kDjDk)  p  (2.2.7) 

for  some  Xk  >  0,  where  the  matrix  Hk  +  XkDjDk  is  at  least  positive  semi-definite. 

In  practice,  it  has  been  found  to  be  more  satisfactory  to  control  the  value  of  6k 
directly  rather  than  A*,  (see  More  [1983]).  Increases  and  decreases  in  6k  are  usually 
based  on  comparing  the  actual  reduction 

?{xk  +  pk)  -  T{xk) 

to  the  reduction  predicted  by  the  model, 

VfIpk  +  \pTkHkPk. 

The  updating  procedure  for  Sk  can  be  as  simple  as  multiplying  the  current  value  by  some 
prescribed  factor,  without  compromising  global  convergence  properties  (see  below).  The 
preferred  strategy  for  decreasing  Sk  is  more  complicated.  An  approximate  minimum  rk 


9 


of  T{xk  +  Tpk)  is  computed  by  safeguarded  polynomial  interpolation  (as  in  linesearch 
methods  —  see  Section  2.2.1),  and  Tfc||Z)fcPfcjj2  is  taken  to  be  the  new  value  of  6k.  It 
may  be  necessary  to  decrease  6k  a  number  of  times  before  a  suitable  reduction  in  T  is 
achieved  and  the  step  to  a  new  iterate  can  be  taken. 

Once  6k  is  assigned  a  value,  it  remains  to  find  pk  when  the  solution  to  (2.2.6) 
is  not  an  unconstrained  minimum.  More  and  Sorensen  [1983]  obtain  A k  in  (2.2.7)  by 
truncating  a  numerical  procedure  for  computing  a  zero  of  the  function 


*(A)  =  ||p*(A)||a 


-6k, 


(2.2.8) 


based  on  the  work  of  Hebden  [1973]  (see  also  Reinsch  [1971]  and  Gay  [1981]).  The 
algorithm  of  Gay  [1983],  implemented  in  the  PORT  Library  [1984],  approximates  p*( A) 
by  a  linear  combination  of  the  (scaled)  steepest  descent  direction  and  the  Newton 
direction.  This  technique  was  devised  by  Powell  [1970]  (see  also  Dennis  and  Mei  [1979]), 
and  is  used  because  it  achieves  comparable  performance  to  methods  that  attempt  to 
approximate  \P(A)  closely,  with  considerably  less  computational  effort. 

Somewhat  stronger  convergence  results  have  been  proven  for  trust-region  meth¬ 
ods  than  are  known  for  linesearch  methods  (see  Section  2.2.1).  Trust-region  methods 
converge  to  an  isolated  local  minimum  under  fairly  mild  conditions  when  exact  second 
derivatives  are  used,  and  otherwise  to  a  stationary  point.  Although  global  convergence 
properties  are  not  affected,  in  practice  the  choice  of  6q  and  the  updating  strategy  for  6k 
are  important.  As  6k,  and  hence  the  norm  of  p,  is  made  to  approach  zero,  the  minimizer 
of  the  quadratic  becomes  parallel  to  the  steepest  descent  direction,  -V.Ffc.  Small  values 
of  6k  are  accordingly  safe,  in  the  sense  that  they  guarantee  a  decrease,  but  progress  may 
be  unacceptably  slow  if  no  provision  is  made  for  taking  larger  steps  where  possible. 

For  more  detail  and  general  discussion  of  trust-region  methods,  see  More  [1983], 
Shultz,  Schnabel,  and  Byrd  [1985],  and  Bulteau  and  Vial  [1987].  Related  variants  are 
described  in  Bulteau  and  Vial  [1985]  and  Byrd,  Schnabel,  and  Shultz  [1988]. 


2.2.3  Stationary  Points  and  Directions  of  Negative  Curvature 
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It  is  possible  to  decrease  T  at  a  stationary  point  x*  if  the  Hessian  matrix  has  one 
or  more  negative  eigenvalues.  The  decrease  is  obtained  by  moving  along  a  direction  of 
negative  curvature;  in  other  words,  a  direction  p  for  which  pTV2Jr(x')p  <  0.  Trust- 
region  methods  that  use  the  quadratic  model  with  exact  Hessian  information  (see  Section 
2.3.1)  will  yield  directions  of  negative  curvature  at  stationary  points  when  V2T(x*)  is 
indefinite,  whereas  the  linesearch  methods  discussed  above  terminate  when  the  gradient 
vanishes. 

A  fundamental  problem  is  that  of  deciding  the  length  of  the  step  to  be  taken 
along  a  direction  of  negative  curvature.  This  problem  is  very  much  related  to  the 
problem  of  setting  a  maximum  step  length  in  order  to  safeguard  a  linesearch  method,  or 
that  of  determining  the  step  bound  in  a  trust-region  method.  First-  and  second-order 
information  about  the  function  at  x*  indicates  that  an  infinite  step  can  be  taken,  since 
the  quadratic  part  of  the  Taylor  series  at  x*  is  unbounded  below  when  V2f(x*)  is 
indefinite.  Clearly  an  infinite  step  is  not  possible  if  T  has  a  finite  minimum. 

Neither  the  question  of  choosing  a  direction  of  negative  curvature,  nor  the  problem 
of  choosing  a  steplength  along  such  directions  has  been  adequately  resolved,  and  thus 
in  most  methods  directions  of  negative  curvature  are  not  explicitly  sought  at  arbitrary 
points.  For  research  on  generating  directions  of  negative  curvature,  and  on  their  use 
in  unconstrained  optimization  algorithms,  see  Gill  and  Murray  [1974a],  Fletcher  and 
Freeman  [1977],  McCormick  [1977],  More  and  Sorensen  [1979],  Goldfarb  [1980],  and 
Shultz,  Schnabel,  and  Byrd  [1985]. 

2.3  Defining  the  Quadratic  Model 

2.3.1  Second-Derivative  Methods 

There  are  two  basic  frameworks  for  defining  Hk  in  the  quadratic  model  (2.1.1)  when 
second  derivative  information  is  available:  direct  modification  of  the  Hessian,  and  trust- 
region  methods.  Both  can  be  viewed  as  procedures  for  producing  a  positive-definite 
quadratic  model  by  modifying  the  exact  Hessian  V2^fc.  A  method  that  combines  the 
two  approaches  is  given  in  Chapter  5  (Section  5)  of  Dennis  and  Schnabel  [1983]. 
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The  modified  Newton  method  of  Giii  and  Murray  [1974a]  is  a  linesearch  method  in 
which  the  definition  of  the  search  direction  is  based  on  the  fact  that  if  Hk  is  positive 
definite,  it  can  be  characterized  by  its  Cholesky  factorization 

Hk  =  RjRk,  (2.3.1) 

where  il*  is  upper-triangular  and  nonsingular  (see,  e.  g.,  Stewart  [1973],  Chapter  3). 
Gill  and  Murray  alter  the  Cholesky  factorization  procedure  so  that  it  can  be  continued 
in  the  event  of  indefiniteness  or  near-singularity.  The  modified  factorization  is  applied  to 
the  Hessian  matrix  V2/*,  resulting  in  the  Cholesky  factorization  of  a  matrix  Hk  with  a 
prescribed  upper  bound  on  its  condition  number.  The  matrix  Hk  may  differ  from  V2Tk 
only  in  the  diagonal  elements.  Local  convergence  properties  of  of  Newton's  method 
are  preserved,  because  Hk  =  V2^  whenever  V2Tk  is  sufficiently  positive  definite. 
An  implementation  is  available  in  the  NAG  Library  [1984]  (subroutine  E04LBF).  For 
information  on  other  direct  modification  methods,  see  Gill,  Murray,  and  Wright  [1981], 
Chapter  4,  and  Higham  [1986]. 

In  trust-region  methods  with  exact  Hessian  information,  a  subproblem  of  the  form 

minVjFjp+|pTV2:Ffcp  (2.3.2) 

subject  to||Z?fcP||2  <  6k, 

is  solved  for  the  step  p*  to  the  next  iterate.  We  recall  from  the  overview  of  trust-regions 
in  Section  2.2.2  that  solving  (2.3.2)  is  equivalent  to  solving 

min  \  pT  (V2;Ffc  +  A kDj Dk)  p  (2.3.3) 

for  some  non-negative  value  of  A*,  with  V2^  +  A kDjDk  positive  semidefinite.  In 
particular,  A*  will  be  positive  whenever  V2^  is  indefinite,  and  also  when  V2^>  is 
positive-definite  if  6k  happens  to  be  smaller  than  the  norm  of  the  scaled  unconstrained 
minimum  of  the  quadratic  objective.  In  contrast  to  the  modified  Newton  method  de¬ 
scribed  above,  all  of  the  eigenvalues  of  V2Tk  are  changed  when  A*  >  0  in  (2.3.3).  As 
long  as  the  constraint  in  (2.3.2)  is  inactive  near  a  local  minimum,  the  local  convergence 
behavior  of  Newton’s  method  is  preserved.  A  recent  implementation  of  a  trust-region 
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method  that  uses  second  derivatives  is  available  in  the  PORT  Library  [1984]  (subroutine 
DMNH;  see  also  Gay  [1983]).  For  further  information  on  trust-regions  with  exact  Hessian 
information,  see  Fletcher  [1980],  Chapter  5,  Gay  [1981],  Sorensen  [1982],  More  [1983], 
More  and  Sorensen  [1983],  and  Shultz,  Schnabel,  and  Byrd  [1985]. 


2.3.2  Quasi-Newton  Methods 


In  quasi-Newton  methods  (also  called  va.ria.ble  metric  or  secant  methods),  a  se¬ 
quence  of  approximations  to  the  Hessian  matrix  of  T  is  generated,  with 

Hk+ i  depending  on  Hk  as  well  as  on  gradient  information  at  the  current  iterate.  The 
approximate  Hessian  is  required  to  satisfy  the  quasi-Newton  condition 

-tffc+i-s*  =  1/fc*  (2.3.4) 


Sk  =  Xk+ 1  -  Xk\  Vk  =  Vfk+l  ~  V/fc. 


The  quantity  yjsk  approximates  the  curvature,  of  T  along  s k.  Equation 

(2.3.4)  does  not  uniquely  define  Hk+\,  papers  that  discuss  completion  of  the  specifica¬ 
tion  include  Dennis  and  More  [1977],  Nazareth  [1984],  T odd  [1984],  and  Flachs  [1986]. 
Conditions  imposed  on  the  approximate  Hessian  almost  always  include  symmetry  and 
positive  definiteness. 

It  is  generally  agreed  that  the  best  overall  performance  is  achieved  by  the  BFGS 
update 


Hk+i  =Hk- 


Hksk{Hksk)T  VkVk 


sjHkSk  yjsk' 

although  precise  reasons  for  its  superiority  are  still  not  known  (see,  e.  g.,  Brodlie  [1977]). 
Like  most  proposed  updates,  the  BFGS  update  is  a  rank-two  modification  of  the  cur¬ 
rent  approximate  Hessian.  The  BFGS  update  preserves  positive  definiteness  whenever 
yjsk  >  0,  a  condition  that  holds  automatically  in  a  linesearch  method  satisfying  (2.2.4). 

Originally,  quasi- Newton  updates  were  formulated  in  terms  of  H^x  rather  than  Hk, 
so  that  minimizing  the  quadratic  (2.1.1)  at  each  stage  in  a  linesearch  algorithm  involved 
only  a  matrix  multiplication  ( 0(n 2)  arithmetic  operations)  rather  than  solution  of  a 
linear  system  (C)(n3)  arithmetic  operations).  Gill  and  Murray  [1972]  showed  that  rank- 
two  quasi-Newton  methods  could  be  implemented  in  0(n2)  operations  per  iteration 
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by  applying  an  update  to  a  Cholesky  factor  (see  Section  2.3.1)  of  Hk-  This  has  the 
additional  advantage  that  it  allows  the  numerical  positive  definiteness  of  Hk  to  be 
monitored  from  iteration  to  iteration.  For  more  information  on  computational  aspects 
of  the  update,  see  Dennis  and  Schnabel  [1983],  Chapter  9,  and  Gill  et  al.  [1985]. 

The  BFGS  method  belongs  to  a  class  of  quasi- Newton  methods  that  can  be  derived 
by  minimizing  the  difference  ( Hk+\  -  Hk)  or  —  *),  'n  various  weighted  norms, 

subject  to  (2.3.4)  [Dennis  and  Schnabel  (1979)].  Other  classes  of  methods  attempt  to 
minimize  the  condition  number  of  Hk  by  selecting  parameters  in  a  class  of  updates  at 
each  step  [Shanno  and  Kettler  (1970);  Oren  (1973,  1982);  Davidon  (1975);  Oren  and 
Spedicato  (1976);  Spedicato  (1976);  Schnabel  (1978)].  Al-Baali  and  Fletcher  [1985] 
apply  a  scaling  factor  before  updating  that  minimizes  an  approximate  measure  of  the 
error  in  the  inverse  Hessian  matrix.  Performance  tests  indicate  that  these  modified 
methods  are  not  as  successful  as  the  BFGS  method  for  general  problems  [Brodlie  (1977); 
Shanno  and  Phua  (1978b);  Al-Baali  and  Fletcher  (1985)]. 

Under  the  same  assumptions  as  required  for  local  quadratic  convergence  of  Newton's 
method,  quasi-Newton  methods  are  locally  superlinearly  convergent,  provided  Ho  is 
sufficiently  close  to  V2Jr(xo)  [Broyden,  Dennis,  and  Mord  (1973)].  Selection  of  the 
initial  Hessian  approximation  Ho  can  be  critical  in  the  performance  of  a  quasi-Newton 
method.  Often  the  identity  is  chosen  because  it  gives  the  steepest-descent  direction 
on  the  first  iteration,  and  it  is  positive  definite.  Computational  tests  have  shown  that 
improved  performance  can  sometimes  be  achieved  by  scaling  Ho  before  performing 
the  first  update  [Shanno  and  Phua  (1978a);  Dennis  and  Schnabel  (1983),  Chapter 
9].  Another  possibility  is  to  use  a  finite-difference  approximation  to  V2.F(zo)  for  Ho, 
modified  if  necessary  to  ensure  positive  definiteness.  Although  the  choice  of  Ho  can 
have  a  significant  effect  on  performance,  the  question  of  how  best  to  choose  Ho  is  still 
open.  It  is  generally  agreed  that  exact  or  approximate  curvature  information  should 
be  used  to  start  the  algorithm  if  it  is  available  at  a  reasonable  cost.  For  a  nonlinear 
least-squares  problem,  Jo  Jo  can  be  used  as  the  initial  estimate,  provided  the  columns 
of  Jo  are  linearly  independent. 
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3.  Methods  for  Nonlinear  Least  Squares 


3.1  Gauss-Newton  Methods 

The  Gauss-Newton  method  is  a  linesearch  method  in  which  the  search  direction  at 
the  current  iterate  minimizes  the  quadratic  function 

9TP+^PTJTJP-  (3.1.1) 

As  a  model  for  the  change  in  the  least-squares  objective,  (3.1.1)  has  the  advantage  that 
it  involves  only  first  derivatives  of  the  residuals,  and  that  JTJ  is  always  at  least  positive 
semi-definite.  If 

paN  =  arg  min  gTp^pTJrJp, 

p£R“  i 

then 

JTJpaN  =  - g ,  (3.1.2) 

so  that  paN  is  a  direction  of  descent  for  fTf  whenever  <7  /  0,  as  required  in  a  linesearch 
method.  To  guarantee  convergence,  the  sequence  of  search  directions  must  also  be 
bounded  away  from  orthogonality  to  the  gradient,  a  condition  that  may  not  be  met  by 
successive  Gauss-Newton  directions  unless  the  eigenvalues  of  JT J  are  bounded  away 
from  zero.  Powell  [1970a]  gives  an  example  of  convergence  of  a  Gauss-Newton  method 
with  exact  linesearch  to  a  non-stationary  point. 

The  Gauss-Newton  method  can  be  viewed  as  a  modification  of  Newton's  method 
in  which  JT  J  is  used  to  approximate  the  Hessian  matrix 

m 

jTj  +  ]C^v2<k  =  jTj  +  b 
1=1 

of  the  nonlinear  least-squares  objective  function.  The  assumption  is  that  the  matrix  JTJ 
should  be  a  good  approximation  to  the  full  Hessian  when  the  residuals  are  small.  In 
fact,  if  f(x*)  =  0  and  J(x*)T J(x*)  is  positive  definite,  then  the  sequence  { x k  +  p%N} 
is  locally  quadratically  convergent  to  z*.  because  J(ik)TJ(xk )  is  an  0{\\xk  -  x*||) 
approximation  to  the  Hessian  of  the  nonlinear  least-squares  objective  at  x*. 
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When  Jt<7  is  singular,  or,  equivalently,  when  J  has  linearly  dependent  columns, 
(3.1.1)  does  not  have  a  unique  minimizer.  The  set  of  vectors  that  minimize  (3.1.1)  is 
the  same  as  the  set  of  solutions  to  the  linear  least-squares  problem 


min  ||  Jp  -f  /||2  ■  (3.1.3) 

One  (theoretically)  well-defined  alternative  that  is  often  approximated  computationally 
is  to  require  the  unique  solution  of  minimum  1%  norm: 


min||p||2,  (3.1.4) 

where  S  is  the  set  of  solutions  to  (3.1.3),  while  another  is  to  replace  J  in  (3.1.3)  by  a 
maximal  linearly  independent  subset  of  its  columns.  In  finite-precision  arithmetic,  there 
is  often  some  ambiguity  about  how  to  formulate  and  solve  an  alternative  to  (3.1.3) 
when  the  columns  of  J  are  "nearly”  linearly  dependent,  which  is  significant  because 
the  numerical  solution  of  these  problems  is  dependent  on  the  criteria  used  to  estimate 
the  rank  of  J.  Fraley  [1987b]  gives  some  detailed  examples  that  illustrate  some  of  the 
difficulties  that  arise  in  implementation. 


3.2  Levenberg-Marquardt  Methods 

Levenberg-Marquardt  methods  alter  the  Gauss-Newton  search  direction  in  the  range 
of  J,  by  replacing  JT  J  in  the  quadratic  model  function  with  JTJ  +  \DTD,  where  A  >  0 
and  D  is  a  diagonal  scaling  matrix  with  positive  diagonal  entries.  The  step  p  between 
successive  iterates  minimizes  the  quadratic  model 

9TP  +  \  PT(JTJ  +  XDTD)p ,  (3.2.1) 

for  some  A  >  0.  Since  the  matrix  JTJ  +  A DT D  is  positive  semidefmite,  minimizers  p\ 
of  (3.2.1)  satisfy  the  equations 

(JrJ  -I-  A DTD)p  =  - g =  -JT/,  (3.2.2) 


which  are  the  normal  equations  for  the  linear  least-squares  problem 

|(n/Ab)','(o)II 


mm 

p6»n 


(3.2.3) 
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Equivalently,  p  solves 


(3-2.4) 


mmgTp+±pTj'IJp 

subject  to  ||Z?p||2  <  6, 

for  some  6  >  0;  that  is,  the  Gauss-Newton  quadratic  model  is  minimized  subject  to  a 
trust-region  constraint. 

Considerable  research  effort  has  been  directed  toward  improvements  in  this  class  of 
methods  since  their  introduction  by  Levenberg  [1944],  Morrison  [1960],  and  Marquardt 
[1963]  for  nonlinear  least-squares  problems,  and  independently  by  Goldfeld,  Quandt,  and 
Trotter  [1966]  for  general  unconstrained  optimization.  More  [1978]  gives  an  implemen¬ 
tation  in  which  he  adjusts  the  step  bound  6  in  (3.2.4)  rather  than  A,  a  strategy  used 
in  trust-region  methods  for  unconstrained  optimization  (see  More  [1983]  for  a  survey). 
Changes  in  6  depend  on  agreement  between  the  actual  reduction  in  the  sum  of  squares 

\  (»/(*+ W)lll-  11/(4)11’). 

with  the  reduction 

9TPx  +  \p\JTJP\ 

predicted  by  the  model  JrJ  +  A DT D,  which  is  the  optimum  value  of  the  objective  in 
(3.2.4).  Increases  are  accomplished  by  taking  6k+i  =  2||Z?jtP*:|j2  ►  while  6  is  decreased 
by  multiplying  by  the  factor  7  <  1.  In  order  to  obtain  A  when  the  bound  in  (3.2.4)  is 
active,  the  nonlinear  equation 

*(A)  =  \\Dpx\\2  -  6  =  ||(Jt7  +  A DJD)~1  s|2  -6  =  0  (3.2.5) 

is  approximately  solved  by  truncating  a  safeguarded  Newton  method  based  on  the  work 
of  Hebden  [1973]  (see  also  Reinsch  [1971]).  More  reports  that,  on  the  average,  (3.2.5) 
is  solved  fewer  than  two  times  per  iteration.  He  also  proves  global  convergence  to  a 
stationary  point  of  /T/,  without  assuming  boundedness  for  {A*}.  Many  computational 
details  are  given,  including  an  efficient  method  for  calculating  the  derivative  of  'P(A) 
in  (3.2.5)  that  uses  the  QR  factorization  of  J.  Equation  (3.2.2)  is  solved  by  a  modi¬ 
fication  of  the  two-stage  factorization  described  by  Osborne  [1972]  that  allows  column 
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pivoting.  Subroutine  LMDER  in  MINPACK  [More,  Garbow,  and  Hillstrom  (1980)]  is  an 
implementation  of  the  method.  Variables  are  scaled  internally  in  LMDER  according  to 
the  following  scheme:  the  initial  scaling  matrix  Do  is  the  square  root  of  the  diagonal  of 
JT  J  evaluated  at  xo,  and  the  ith  diagonal  element  of  Dk  is  taken  to  be  the  maximum 
of  the  ith  diagonal  element  of  Dk- 1  and  the  square  root  of  the  tth  diagonal  element 
of  JrJ.  Numerical  results  are  presented  indicating  that  this  scaling  compares  favorably 
with  those  used  in  earlier  research.  The  user  also  has  the  option  of  providing  an  initial 
diagonal  scaling  matrix  that  is  retained  throughout  the  computation. 

3.3  Corrected  Gauss-Newton  Methods 

Gill  and  Murray  [1976]  propose  a  linesearch  algorithm  that  divides  into  comple¬ 
mentary  subspaces  7Z  and  AT,  where  R  C  R(  JT),  and  J\f  is  nearly  orthogonal  to  R(JT). 
The  search  direction  is  the  sum  of  a  Gauss- Newton  direction  in  R,  and  a  projected  New¬ 
ton  direction  in  A T.  This  strategy  avoids  a  shortcoming  of  Gauss- Newton  methods  — 
that  components  of  the  search  direction  that  are  nearly  orthogonal  to  R(JT)  may  not 
be  well  determined  when  J  is  ill-conditioned  —  because  each  component  is  computed 
from  a  reasonably  well-conditioned  subproblem.  The  vector  x  -  x*  may  become  almost 
entirely  in  R(JT)  in  a  Gauss-Newton  method,  yet  the  algorithm  computes  a  search 
direction  that  is  virtually  orthogonal  to  R(JT)  due  to  ill  conditioning  in  the  Jacobian 
(see  Fraley  [1987b]).  Gill  and  Murray  show  that  both  Gauss-Newton  algorithms  defined 
by  (3.1.4)  and  Levenberg-Marquardt  algorithms  generate  search  directions  that  lie  in 
R(Jr),  while  the  Newton  search  direction  generally  will  have  a  component  in  A f{J), 
the  orthogonal  complement  of  R(JT),  whenever  J  has  linearly  dependent  columns.  For 
problems  with  small  residuals,  they  point  out  that  JTJ  is  a  reasonable  approximation 
to  the  full  Hessian  in  7 Z(JT),  but  not  in  Af(J).  Thus,  in  situations  where  x  -  x*  is 
orthogonal  to  R(Jr),  and  J  is  well-conditioned  but  has  linearly  dependent  columns  (for 
example,  when  m  <  n),  the  Gauss-Newton  and  Levenberg-Marquardt  directions  have 
no  component  in  the  direction  of  x  —  z*,  while  Newton’s  method  and  also  the  method 
of  Gill  and  Murray  would  have  components  in  both  R(JT)  and  M{J). 
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A  version  of  this  algorithm  called  the  corrected  Gauss-Newton  method  [Gill  and 
Murray  (1978)]  forms  the  basis  for  the  nonlinear  least-squares  software  in  the  NAG 
Library  [1984],  Rules  based  on  the  relative  size  of  the  singular  values  of  J  are  given 
for  choosing  an  integer  grade(J)  to  approximate  rank(J),  and  an  attempt  is  made 
to  group  together  singular  values  that  are  similar  in  magnitude.  The  method  is  not 
as  sensitive  to  grade(J)  as  Gauss-Newton  is  to  rank  estimation,  both  because  of  the 
division  of  the  computation  of  the  search  direction  into  separate  components  in  7Z  and 
M,  and  because  grade(J)  is  varied  adaptively  based  on  a  measure  of  the  progress  of 
the  minimization.  Moreover,  the  rate  of  convergence  is  potentially  faster  than  Gauss- 
Newton  or  Levenberg-Marquardt  methods  on  problems  with  nonzero  residuals.  The 
quantity  grade(J)  is  reduced  when  the  sum  of  squares  is  not  adequately  decreasing,  so 
that  there  is  the  potential  of  having  M  =  3?n  (with  exact  second  derivatives,  this  implies 
taking  full  Newton  steps)  in  the  vicinity  of  a  solution. 

When  rank(J)  =  grade(J)  =  n,  the  search  direction  p  is  a  full-rank  Gauss- 
Newton  direction.  Otherwise  the  vector  p  is  computed  as  the  sum  of  two  mutually 
orthogonal  components:  a  Gauss- Newton  direction,  and  a  projected  Newton  direction. 
The  projected  Hessian  is  replaced  by  a  modified  Cholesky  factorization  (see  Section 
2.3.1)  if  it  is  computationally  singular  or  indefinite.  A  modified  Newton  search  direction 
(corresponding  to  the  case  grade(J)  =  0)  is  used  whenever  if  |cos(<7,p)|  is  smaller  than 
some  prescribed  value,  or  if  gTp  is  positive.  A  quasi-Newton  approximation  to  B  (see  the 
discussion  in  Section  3.4)  and  a  finite-difference  approximation  to  the  projected  matrix 
ZT  BZ  along  the  columns  of  Z,  where  Z  is  an  orthogonal  basis  for  M,  are  given  as 
alternatives  to  handle  cases  in  which  second  derivatives  of  the  residual  functions  are  not 
available  or  are  difficult  to  compute.  Gill  and  Murray  test  their  method  on  a  set  of  twenty- 
three  problems,  and  find  that  the  version  that  uses  quasi-Newton  approximations  to  B 
does  not  perform  as  well  as  those  that  use  exact  second  derivatives  or  finite-difference 
approximations  to  a  projection  of  B.  They  observe  only  linear  convergence  for  the  quasi- 
Newton  version  on  problems  with  large  residuals.  The  algorithms  are  implemented  in  the 
NAG  Library  [1984];  subroutine  E04HEF  uses  exact  second  derivatives,  while  subroutine 
E04GBF  is  the  quasi-Newton  version. 
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3.4  Special  Quasi-Newton  Methods 

Special  quasi-Newton  methods  for  nonlinear  least  squares  use  a  Hessian  of  the  form 
JT  J  -|-  B  in  the  quadratic  model,  so  that  the  search  direction  differs  from  the  Gauss- 
Newton  direction  in  TZ{  JT),  and  also  has  a  component  in  Af(  J)  when  J  is  rank-deficient. 
The  matrix  B  is  a  quasi-Newton  approximation  to  the  term  B  =  4>C^24>i  in  the 

Hessian  of  the  nonlinear  least-squares  objective.  Brown  and  Dennis  [1971]  first  proposed 
a  method  in  which  the  Hessian  matrix  of  each  of  the  residuals  was  updated  separately. 
This  approach  is  impractical  because  it  entails  the  storage  of  m  symmetric  matrices  of 
order  n,  and  more  recent  research  has  aimed  to  approximate  B  as  a  sum. 

Gill  and  Murray  [1978]  discuss  a  linesearch  method  in  which  they  use  the  augmented 
Gauss-Newton  quadratic  model  only  to  compute  a  component  of  the  search  direction 
in  a  subspace  that  approximates  the  null  space  of  the  Jacobian  (see  Section  3.3).  They 
apply  the  BFGS  formula  for  unconstrained  optimization  (see,  e.  g.,  Dennis  and  More 
[1977])  to  the  matrix  Hk  =  «/J+1  Jk+i  +  Bk  with  the  quasi-Newton  condition 

Hk+lSk  =  Vk 

where 

Sk  =  z/t+i  -  xfc  and  y*  =  gk+i  -  9k , 

and  then  form  Bk+i  =  H k+i  -  Jj+i^k+i-  They  point  out  that,  if  Jj+i^k+i  +  's 
positive  definite,  and  yjs*  >  0,  then  Jj+i^k+i  +  Bk+i  is  also  positive  definite  with 
this  scheme.  In  order  to  safeguard  the  method,  the  projected  approximate  Hessian 
is  replaced  by  a  modified  Cholesky  factorization  when  it  is  singular  or  indefinite.  In 
addition,  if  the  cosine  of  the  angle  between  the  search  direction  and  the  gradient  of  the 
nonlinear  least-squares  objective  exceeds  a  fixed  threshold  value,  a  modified  Newton  step 
with  the  full  augmented  approximate  Hessian  is  taken.  See  Section  3.3  for  a  summary 
of  their  observations  on  the  performance  of  the  methods. 

Dennis,  Gay,  and  Welsch  [1981a]  apply  a  scaled  DFP  update  to  Bk  at  each  step. 
The  new  approximation  Bk+\  solves 
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subject  to 


Hsic  =  j ik]  H  positive  definite 

Bsk  =  Jj+1fk+i  -  Jkfk+i;  B  symmetric, 


where 

r*  = 

The  scale  factor  r*  is  based  on  the  observation  that  the  quasi- Newton  approximation 
to  B  is  often  too  large  with  the  unsealed  update,  on  account  of  the  contribution  of 
the  residuals.  The  term  yJsk/sjBkSk  in  r*  is  derived  from  the  self-scaling  principles 
for  quasi-Newton  methods  of  Oren  [1973],  and  attempts  to  shift  the  eigenvalues  of  the 
approximation  Bk  to  overlap  with  those  of  Bk,  using  curvature  information  at  i* .  The 
algorithm  forms  the  basis  for  the  ACM  computer  program  NL2S0L  [Dennis,  Gay,  and 
Welsch  (1981b)],  which  is  distributed  by  the  PORT  Library  [1984]  as  subroutines  N2G 
and  DN2G.  It  is  implemented  as  an  adaptive  method,  in  that  Gauss-Newton  steps  are 
taken  if  the  Gauss-Newton  quadratic  model  predicts  the  reduction  in  the  function  better 
than  the  quadratic  model  that  includes  the  term  involving  B.  A  trust-region  strategy 
is  used  to  enforce  global  convergence.  Numerical  results  are  given  in  Dennis,  Gay,  and 
Welsch  [1981a]  for  a  set  of  twenty-four  test  problems,  many  with  two  or  three  different 
starting  values. 
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4.  Discussion  and  Summary  of  Numerical  Results 


In  this  section,  we  summarize  numerical  results  obtained  for  the  unconstrained 
optimization  and  nonlinear  least-squares  methods;  more  detailed  results  are  tabulated 
in  the  Appendix.  The  tests  were  performed  using  the  following  software: 

•  DMNG/SUMSOL  -  Trust-region  method  for  unconstrained  optimization  that  uses  a 
quasi- Newton  approximation  to  the  Hessian  matrix.  From  the  PORT  Library  [1984], 

•  NPSOL  -  Linesearch  method  for  unconstrained  optimization  that  uses  a  quasi-Newton 
approximation  to  the  Hessian  matrix.  From  the  Systems  Optimization  Laboratory, 
Stanford  University  (see  Gill  et  al.  [1986b;  1987]).  Also  available  in  the  NAG 
Library. 

•  DMNH/HUMSOL  -  Trust-region  method  for  unconstrained  optimization  that  uses  ana¬ 
lytic  second  derivatives.  From  the  PORT  Library  [1984], 

•  MNA  -  Linesearch  method  for  unconstrained  optimization  that  uses  analytic  second 
derivatives.  This  implementation,  which  is  available  at  Stanford  Linear  Accelerator 
Center,  is  from  the  National  Physical  Laboratory,  England.  It  is  essentially  the  same 
as  subroutine  E04LBF  from  the  NAG  Library  [1984], 

•  G-N  -  Gauss-Newton  algorithm  for  nonlinear  least  squares  that  uses  LSSOL  (Gill  et 
al.  [1986a])  to  solve  the  linear  least  squares  subproblems,  and  the  linesearch  from 
NPSOL  (Gill  et  al.  [1986b]).  Both  LSSOL  and  NPSOL  are  also  available  in  the  NAG 
Library. 

•  LMDER-  Levenberg-Marquardt  method  for  nonlinear  least  squares.  From  MINPACK 
(More,  Garbow,  and  Hillstrom  [1980]). 

•  DN2G/NL2S0L  -  Adaptive  method  for  nonlinear  least  squares  (combines  Gauss- 
Newton,  Levenberg-Marquardt,  and  special  quasi-Newton  methods).  From  the 
PORT  Library  [1984], 

•  LSQFDQ  -  Corrected  Gauss-Newton  method  that  uses  a  quasi-Newton  approxima¬ 
tion  to  the  Hessian  matrix.  This  implementation,  which  is  available  at  Stanford 
Linear  Accelerator  Center,  is  from  the  National  Physical  Laboratory,  England.  It  is 
essentially  the  same  as  subroutine  E04GBF  from  the  NAG  Library  [1984], 
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•  LSQSDN  -  Corrected  Gauss-Newton  method  that  uses  analytic  second  derivatives. 
This  implementation,  which  is  available  at  Stanford  Linear  Accelerator  Center,  is 
from  the  National  Physical  Laboratory,  England.  It  is  essentially  the  same  as  sub¬ 
routine  E04HEF  from  the  NAG  Library  [1984]. 

Information  about  the  individual  test  problems  is  given  in  the  Appendix.  The  number 
of  function  evaluations  required  by  each  subroutine  is  listed  in  the  tables  at  the  end  of 
this  section.  In  addition,  the  following  symbols  are  used: 

0  -  zero-residual  problem 

1  -  linear  least-squares  problem 

-  -  failure  to  achieve  an  approximate  solution 

appears  to  be  unable  to  terminate  at  an  approximate  solution 
1  -  local  minimum 

*  -  termination  criteria  satisfied  at  a  point  away  from  a  local  minimum 

*  -  failed  with  default  step  length  or  trust-region  size 

Two  columns  of  figures  corresponding  to  two  different  values  of  a  single  parameter 
are  given  for  each  subroutine.  For  the  Gauss-Newton  methods,  the  parameter  affects 
rank  estimation;  for  all  of  the  other  methods,  the  parameter  affects  termination  criteria. 
See  the  tables  of  numerical  results  given  in  the  Appendix  for  information  about  the 
precise  choices  that  were  made.  The  wide  variability  in  the  numerical  results  makes  it 
difficult  to  draw  definitive  conclusions  about  the  relative  performance  of  the  software, 
because  observations  of  small  samples  could  result  in  misleading  generalizations.  The 
sources  of  this  variability  are  discussed  below. 

First,  the  number  of  function  evaluations  may  not  be  an  adequate  basis  for  compar¬ 
ison.  The  routines  vary  in  the  number  of  gradient  evaluations  performed  per  function 
evaluation,  and  second-derivative  methods  require  evaluation  of  the  Hessian  matrix. 
Moreover,  when  function  evaluations  are  relatively  inexpensive,  costs  could  be  domi¬ 
nated  by  other  portions  of  the  computation.  Another  difficulty  in  making  comparisons 
is  that  the  definition  of  an  acceptable  minimum  varies  from  routine  to  routine.  For 
example,  the  norm  of  the  gradient  of  the  nonlinear  least-squares  objective,  ||$||,  at  an 
alleged  solution  x*  may  differ  considerably  for  different  software,  although  <?(x*)  =  0  is 
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a  necessary  condition  for  a  minimum.  (On  problem  10.,  LMDER  terminates  at  a  point 
for  which  j|^||  is  of  order  10°,  while  DN2G  terminates  at  a  point  for  which  ||^||  is  of  order 
10~3.)  Most  algorithms  do  not  attempt  to  reduce  ||p||  directly,  but  convergence  criteria 
may  include  a  threshold  on  ||<7(a:*)||.  Depending  on  how  this  threshold  is  used  in  relation 
to  other  criteria,  some  routines  may  spend  more  function  evaluations  in  anticipation  of 
a  reduction  in  ||p||  than  others.  A  small  value  of  ||g||  means  greater  certainty  that  a 
minimum  has  actually  been  obtained,  but  may  be  unreasonably  expensive  to  achieve  in 
practice. 

Second,  aside  from  design  choices  that  define  a  particular  implementation  of  an 
algorithm,  the  user  is  permitted  to  specify  certain  parameters  that  may  affect  perfor¬ 
mance.  Fraley  [1987b]  gives  examples  that  illustrate  the  sensitivity  of  Gauss-Newton 
methods  to  rank  estimation  criteria  (see,  e.  g.,  problems  35b.,  36a.,  and  20d.).  For 
problems  on  which  an  algorithm  is  linearly  convergent,  small  changes  in  tolerances  that 
are  used  to  define  convergence  criteria  can  mean  substantial  differences  in  the  amount 
of  computation  required  in  order  to  obtain  a  point  satisfying  conditions  for  convergence 
(see,  e.  g.,  DMNG  on  24b.,  and  LMDER  on  40.).  Selection  of  a  maximum  steplength  or 
an  initial  trust-region  radius  can  also  be  critical  factor  in  the  performance  of  a  method. 
In  these  tests,  the  default  values  for  these  parameters  were  altered  only  in  cases  where  a 
method  was  initially  observed  to  fail  by  attempting  to  evaluate  problem  functions  outside 
the  region  in  which  they  are  numerically  defined  (see,  e.  g.,  the  results  for  the  DeVilliers 
and  Glasser  test  problems  42.  and  43.).  Failures  of  this  sort  may  be  caused  either  by 
poor  scaling  among  the  variables  in  the  problem,  or  by  ill-conditioning  within  subprob- 
lems.  Hiebert  [1981]  is  of  the  opinion  that  failure  due  to  ill-conditioning  can  be  avoided 
in  software,  but  that  it  is  not  always  possible  to  anticipate  abnormal  terminations  that 
are  caused  by  bad  scaling. 

In  NPSOL,  one  can  specify  bounds  on  the  variables,  as  well  as  adjust  the  maximum 
step  length,  in  order  to  deal  with  this  type  of  difficulty.  Subroutine  E04LBF,  the  version 
of  MNA  that  is  available  in  the  NAG  Library  [1984],  also  provides  for  bounds  on  variables, 
and  there  are  alternative  versions  of  all  of  the  PORT  software  used  in  these  tests  that 
allow  bounds  to  be  specified.  Unfortunately,  when  bounds  on  the  variables  are  included 
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in  the  problem  formulation,  local  minima  at  which  the  bounds  are  active  may  be  found 
rather  than  local  minima  for  the  nonlinear  least-squares  problem.  See  Hanson  and  Krogh 
[1987]  for  numerical  results  in  which  simple  bounds  are  included  for  some  problems.  Holt 
and  Fletcher  [1979]  give  an  algorithm  designed  for  nonlinear  least-squares  problems  with 
explicit  bounds  on  the  variables. 

Third,  the  performance  of  any  given  method  over  the  set  of  test  problems  is  by  no 
means  uniform,  and  it  is  not  easy  to  separate  the  problems  into  classes  for  which  the 
behavior  of  an  algorithm  can  be  categorized.  One  reason  for  this  is  that  many  of  the  test 
problems  recur  in  the  literature  precisely  because  they  have  certain  distinguishing  prop¬ 
erties.  Powell's  singular  function  and  variants  (13.  and  22.)  are  zero-residual  problems 
in  which  the  Jacobian  becomes  singular  at  the  solution.  The  McKeown  test  problems 
(39.,  40.,  and  41.)  are  chosen  so  that  the  Jacobian  is  well-conditioned  everywhere, 
and  the  rate  of  convergence  for  the  Gauss-Newton  method  with  unit  steplength  can  be 
controlled  by  varying  a  single  parameter  (the  parameter  can  also  be  chosen  so  that  the 
unit-step  Gauss-Newton  method  diverges).  Both  Powell’s  singular  function  and  McK- 
eown's  test  problems  are  constructed  analytically  rather  than  derived  from  data-fitting 
applications.  The  matrix  square  root  problems  (36.)  are  examples  of  small,  dense, 
nonlinear  systems  of  equations  requiring  a  very  accurate  solution.  Watson's  function 
(20.)  comes  from  polynomial  interpolation,  and  has  multiple  local  minima  with  small, 
but  nonzero,  residuals.  It  also  has  the  feature  that  the  Jacobian  becomes  increasingly 
ill-conditioned  as  the  problem  size  is  increased  (see  Fraley  [1987b]).  The  Gulf  Research 
and  Development  function  (11.)  has  discontinuities  in  the  derivative  of  each  residual 
on  a  one-dimensional  subspace  and  hence  violates  the  assumption  (made  in  developing 
all  of  the  algorithms  we  have  discussed)  that  the  sum  of  squares  has  continuous  second 
derivatives.  The  results  for  the  DeVilliers  and  Glasser  test  problems  (42.  and  43.) 
illustrate  variability  in  performance  due  to  the  use  of  different  starting  values.  More 
generally,  the  behavior  of  a  given  method  for  a  certain  type  of  residual  function  may 
not  be  uniform  over  several  sets  of  defining  data  of  similar  magnitude,  as  shown  by  the 
results  for  the  Dennis,  Gay,  and  Vu  test  problems  (44.  and  45.). 
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Finally,  there  is  considerable  diversity  in  performance  among  the  routines  tested, 
and  few  generalizations  are  possible.  Our  data  generally  supports  the  use  of  nonlinear 
least-squares  software  over  that  designed  for  general  unconstrained  minimization,  but 
there  are  some  exceptions  (see,  e.  g.,  the  McKeown  test  problems  39.  -  41.).  Of  the 
nonlinear  least-squares  routines,  DN2G  (NL2S0L)  is  often  the  best  (the  Dennis,  Gay,  and 
Vu  test  problems  44.  and  45.  are  examples  of  exceptions).  When  second  derivatives 
are  relatively  cheap  to  obtain,  the  use  of  an  unconstrained  optimization  method  that 
uses  exact  second  derivatives  may  be  a  reasonable  alternative  to  a  nonlinear  least- 
squares  method  (see,  e.  g.,  the  results  for  the  penalty  function  23b.).  Our  tests  do 
not  indicate  overall  superiority  of  any  particular  method  over  the  others;  in  situations 
in  which  a  variety  of  problems  are  being  solved,  we  conclude  that  it  is  desirable  to  have 
the  ffexibifty  to  choose  from  among  several  methods. 
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Summary  of  Results  :  Nonlinear  Least-Squares  Methods 


(number  of  function  evaluations) 

Mor£,  Garbow,  and  Hillstrom  Test  Problems 


n  m  DMNG  MPSOL  DMNH  MNA 


l.° 

2 

2 

40 

42 

27 

29 

32 

32 

14 

14 

2.° 

2 

2 

12' 

12' 

9 

10 

10' 

ltf 

8 

8 

3.° 

2 

2 

217 

220 

897 

897 

130 

132 

175 

175 

4.° 

2 

3 

66 

67 

20 

20 

22 

23 

1 

1 

5.° 

2 

3 

16 

17 

20 

22 

11 

12 

35* 

35* 

6. 

2 

10 

33 

34 

14* 

15* 

11 

11 

12 

12 

7.° 

3 

3 

28 

30 

37 

38 

16 

17 

14 

14 

8. 

3 

15 

19 

22 

22 

23 

9 

10 

11 

11 

9. 

3 

15 

8 

12 

8 

9 

4 

4 

3 

3 

10. 

3 

16 

465 

467 

450 

450 

382 

388 

249 

249 

11.° 

3 

4' 

327 

2 ' 

2' 

290 

292 

538 

538 

12.° 

3 

10 

43 

45 

34 

35 

24 

24 

43 

43 

13.° 

4 

4 

62* 

89 

66 

71 

27* 

38 

38 

23*  23* 

14.° 

4 

6 

100 

102 

50 

51 

42 

49 

54 

54 

15. 

4 

11 

35 

36 

33 

35 

11 

12 

20 

20 

16. 

4 

20 

46 

47 

24 

25 

11 

13 

9 

10 

17. 

5 

33 

69 

72 

32*' 

56* 

46* 

47* 

43 

43 

18.° 

6 

13 

45 

47 

45' 

45' 

- 

- 

44 

44 

19. 

11 

65 

69 

72 

88 

90 

23 

24 

7* 

8* 

20a. 

6 

31 

35 

37 

43 

46 

15 

15 

13 

13 

20b. 

9 

31 

76 

79 

83 

85 

20 

22 

14 

14 

20c. 

12 

31 

89' 

148' 

55' 

151 

24 

24 

14 

14 

20d. 

20 

31 

110' 

134' 

73' 

114' 

50' 

*■■*/ 

1295' 

1295' 

2  la.0 

10 

10 

120 

125 

101 

104 

25 

26 

14 

14 

21b.° 

20 

20 

189 

193 

252 

265 

27 

27 

14 

14 

22a.° 

12 

12 

143' 

235 

83* 

165* 

28* 

40 

23* 

23* 

22b.° 

20 

20 

187* 

344 

103* 

196* 

29* 

40 

24* 

24* 

23a. 

4 

5 

77 

78 

198 

198 

42 

43 

43 

43 

23b. 

10 

11 

80 

81 

117 

124 

44 

45 

44 

44 

24a. 

4 

8 

364 

472 

23' 

462 

126 

128 

158 

158 

24b. 

10 

20 

475 

632 

368 

419 

158 

162 

133 

133 
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Summary  of  Results  :  Unconstrained  Optimization  Methods 
(number  of  function  evaluations) 


More,  Garbow,  and  Hillstrom  Test  Problems 

G-N  LMDER  DN2G  LSQFDQ  LSQSDN 


l.° 

31 

31 

22 

22 

14 

14 

31 

31 

31 

31 

2.° 

180' 

235' 

14' 

21' 

10' 

12* 

36' 

36' 

18' 

18' 

3.° 

31' 

42 

19 

19 

64 

65 

47 

47 

47 

47 

4.° 

54 

54 

40' 

54 

40' 

53 

64 

64 

53 

53 

5.° 

8 

8 

9 

10 

9 

10 

14 

14 

10 

10 

6. 

_  1 

_*t 

21 

28 

14 

16 

54 

54 

36* 

36* 

7.° 

13 

13 

11 

12 

13 

14 

20 

20 

14 

14 

8. 

7 

7 

6 

7 

7 

8 

13 

13 

6 

6 

9. 

3 

3 

4 

5 

3 

5 

3 

3 

3 

3 

10. 

30 

126 

126 

132 

133 

18 

18 

17 

17 

11.° 

- 

- 

- 

- 

- 

- 

69' 

69' 

30*' 

30*' 

12.° 

7 

7 

7 

8 

8 

9 

12 

12 

8 

12 

13.° 

16' 

16' 

65 

65 

19' 

25 

18' 

18' 

18' 

18' 

14.° 

96 

96 

70 

70 

52 

53 

81 

81 

87 

99 

15. 

43 

43 

18 

28 

11 

12 

30 

30 

16 

16 

16. 

3651 

3651 

264 

356 

21 

22 

62 

62 

45 

45 

17. 

13 

13 

18 

19 

26 

27 

19 

19 

14 

18 

18.° 

_ f 

- 

46 

46 

45 

46 

_ 9 

_ 9 

243* 

247* 

19. 

24 

24 

17 

19 

20 

22 

33 

33 

19 

19 

20a. 

12 

12 

8 

10 

13 

13 

32 

32 

9 

9 

20b. 

6 

6 

9 

10 

12 

15 

6 

6 

6 

6 

20c. 

6 

6 

10 

12 

14 

14 

6 

6 

6 

6 

20d. 

6' 

6 

18 

23 

7' 

18 

18 

10 

12 

21a.° 

31 

31 

22 

22 

27 

27 

31 

31 

31 

31 

21b.° 

31 

31 

22 

22 

16 

16 

31 

31 

31 

31 

22a.° 

16' 

16' 

72 

72 

20' 

26 

18' 

18' 

18' 

18' 

22b.° 

16' 

16' 

69 

69 

19' 

26 

18' 

18' 

18' 

18' 

23a. 

90 

90 

34 

44 

36 

37 

80 

80 

58 

58 

23b. 

274 

274 

84 

104 

61 

68 

143 

143 

124 

124 

24a. 

1043 

1043 

151 

156 

139 

142 

411 

411 

228 

228 

24b. 

- 

- 

80 

88 

129 

138 

566* 

566* 

150 

150 

28 


Summary  of  Results  :  Nonlinear  Least-Squares  Methods 

(number  of  function  evaluations) 


Mor£,  Garbow,  and  Hillstrom  Test  Problems  (continued) 


n 

m 

DMNG 

NPS0L 

DMNH 

MNA 

25a.° 

10 

12 

20 

21 

19 

20 

15 

15 

14 

14 

25b.° 

20 

22 

25 

26 

24 

24 

18 

19 

17 

18 

26a.° 

10 

10 

34' 

37' 

33' 

35' 

11' 

12' 

21 

21 

26b.° 

20 

20 

62* 

65' 

76' 

78' 

20' 

20' 

30 

30 

27a.° 

10 

10 

13 

16 

18 

19 

9 

10 

22 

22 

27b.° 

20 

20 

15 

18 

30 

32 

11 

12 

31 

31 

28a.° 

10 

10 

31 

34 

33 

36 

4 

4 

4 

4 

28b.° 

20 

20 

60 

64 

54 

56 

4 

4 

4 

4 

29a.° 

10 

10 

8 

10 

7 

8 

4 

5 

4 

4 

29b. 0 

20 

20 

8 

10 

7 

8 

4 

5 

4 

4 

30a.° 

10 

10 

51 

57 

37 

40 

6 

7 

7 

7 

30b.° 

20 

20 

65 

88 

65' 

67' 

6 

7 

7 

7 

31a.° 

10 

10 

46 

60 

67' 

70* 

9 

9 

9 

9 

31b.° 

20 

20 

47 

63 

141 

144 

9 

9 

9 

9 

32.t 

10 

20 

6 

6 

2 

2 

6 

6 

4 

4 

33. 1 

10 

20 

4 

4 

4 

4 

5 

5 

27 

27 

34.1 

10 

20 

5 

5 

4 

4 

6 

6 

20 

20 

35a. 

8 

8 

34 

38 

31 

33 

14 

14 

41 

41 

35b.° 

9 

9 

44 

46 

29 

32 

17 

18 

66 

66 

35c. 

10 

10 

41' 

45' 

37' 

42' 

19 

20 

86 

86 

Matrix  Square  Root  Test  Problems 

n 

m 

DMNG 

NPS0L 

DMNH 

MNA 

36a.° 

4 

4 

- 

- 

_t 

- 

- 

- 

- 

36b.° 

9 

9 

- 

- 

_t 

_1 

- 

- 

- 

- 

36c.° 

9 

9 

69 

101 

3 

3 

31 

35 

3188  3188 

36d.° 

9 

9 

- 

- 

- 

- 

- 

- 

29 


Summary  of  Results  :  Unconstrained  Optimization  Methods 

(number  of  function  evaluations) 


More,  Garbow,  and  Hillstrom  Test  Problems  (continued) 


G-N  LMDER  DN2G  LSQFDQ  LSQSDN 


25a.° 

11 

11 

11 

12 

15 

16 

16 

16 

12 

17 

25b. 0 

13 

13 

13 

14 

19 

19 

18 

18 

14 

19 

26a. 0 

16 

16 

28' 

37' 

11 

12 

22 

22 

18 

22 

26b.° 

20 

20 

57' 

71' 

39 

42 

18 

24 

18 

26 

27a.° 

21 

21 

15 

15 

8 

9 

26’ 

26* 

22* 

28* 

27b.° 

31* 

31* 

5' 

18 

11 

12 

31* 

31* 

21* 

27* 

28a. 0 

4 

4 

5 

5 

4 

4 

4 

4 

4 

4 

28b.° 

4 

4 

5 

5 

3 

4 

4 

4 

4 

4 

29a. 0 

4 

4 

5 

5 

4 

4 

10 

10 

6 

6 

29b.° 

4 

4 

5 

5 

4 

4 

10 

10 

6 

6 

30a.° 

6 

6 

6 

7 

8 

9 

11 

11 

7 

11 

30b. 0 

6 

6 

6 

7 

8 

9 

11 

11 

7 

11 

31a.° 

7 

7 

7 

8 

10 

11 

12 

12 

8 

12 

31b. 0 

7 

7 

7 

8 

10 

11 

12 

12 

8 

12 

32.*- 

2 

2 

3 

3 

5 

5 

2 

2 

2 

2 

33. L 

3 

- 

3 

8 

18 

18 

2 

2 

2 

2 

34. 1 

3 

3 

3 

7 

13 

13 

2 

2 

2 

2 

35a. 

_t 

- 

40 

53 

23 

24 

87* 

87* 

74 

74 

35b. 0 

148 

12 

13 

11 

11 

34 

34 

30 

34 

35c. 

_ 

_f 

25 

34 

17' 

19' 

73*' 

73*' 

43' 

43' 

Matrix  Square  Root  Test  Problems 
G-N  LMDER  DN2G  LSQFDQ  LSQSDN 


36a.° 

2885' 

36 

- 

- 

- 

- 

44 

44 

38 

38 

36b  ® 

683' 

36 

9 

10 

16' 

- 

44 

44 

38 

38 

36c.° 

20 

20 

29 

40 

16 

22 

28 

28 

28 

28 

36d.° 

74' 

- 

2 

2 

4 

4 

424 

424 

380 

380 

30 


Summary  of  Results  :  Nonlinear  Least-Squares  Methods 

(number  of  function  evaluations) 


Hanson  Test  Problems 

n 

m 

DMNG 

NPSOL 

DMNH 

MNA 

37. 

2 

16 

22  22 

14  15 

16  17 

6  6 

38. 

3 

16 

31  32 

21*  23* 

14  14 

13  13 

McKeown  Test  Problems 
DMNG  NPSOL  DMNH  MNA 


39a. 

2 

3 

9 

11 

10 

11 

4 

4 

4 

4 

39b. 

2 

3 

9 

10 

9 

10 

4 

4 

4 

4 

39c. 

3 

6 

7 

6 

8 

4 

5 

4 

4 

39d. 

2 

3 

8 

9 

10 

11 

6 

6 

6 

6 

39e. 

2 

3 

11 

12 

16 

17 

8 

8 

8 

8 

39f. 

2 

3 

11 

11 

28 

29 

11 

11 

11 

11 

39g. 

2 

3 

17 

18 

30 

31 

13 

14 

14 

14 

40a. 

3 

4 

11 

12 

11 

12 

4 

4 

4 

4 

40b. 

3 

4 

10 

12 

11 

12 

4 

5 

4 

4 

40c. 

3 

4 

9 

10 

9 

10 

5 

5 

5 

5 

40d. 

3 

4 

10 

11 

14 

15 

5 

6 

6 

6 

40e. 

3 

4 

11 

13 

19 

20 

7 

7 

8 

8 

40f. 

3 

4 

14 

16 

33 

34 

10 

10 

10 

10 

40g. 

3 

4 

18 

20 

45 

46 

13 

13 

13 

13 

41a. 

5 

10 

11 

13 

12 

12 

4 

4 

4 

4 

41b. 

5 

10 

11 

13 

12 

13 

4 

4 

4 

4 

41c. 

5 

10 

13 

14 

12 

14 

8 

8 

8 

8 

41d. 

5 

10 

17 

20 

17 

20 

8 

9 

9 

9 

41e. 

5 

10 

24 

26 

51 

54 

11 

12 

12 

12 

41f. 

5 

10 

27 

31 

51 

53 

14 

14 

14 

14 

41g. 

5 

10 

32 

35 

62 

69 

17 

17 

17 

17 

31 


Summary  of  Results  :  Unconstrained  Optimization  Methods 

(number  of  function  evaluations) 


Hanson  Test  Problems 

G-N  LMDER  DN2G  LSQFDQ  LSQSDN 


37. 

39 

39 

15 

21 

10 

11 

25 

25 

9 

9 

38. 

58 

58 

18 

28 

10 

12 

30 

30 

10 

10 

McKeown  Test  Problems 

G-N  LMDER  DN2G  LSQFDQ  LSQSDN 


39a. 

8 

8 

5 

6 

5 

5 

17 

17 

4 

4 

39b. 

32 

32 

14 

21 

6 

7 

24 

24 

6 

6 

39c. 

23 

23 

18 

25 

7 

8 

22 

23 

9 

9 

39d. 

681 

681 

20 

28 

7 

8 

31 

32 

12 

12 

39e. 

- 

- 

28 

44 

9 

10 

32 

32 

12 

12 

39f. 

- 

- 

31 

44 

14 

15 

43 

43 

25 

25 

39g. 

- 

- 

39 

44 

18 

20 

49 

49 

39 

39 

40a. 

13 

13 

6 

9 

7 

7 

18 

18 

5 

5 

40b. 

16 

16 

14 

17 

7 

11 

19 

19 

6 

6 

40c. 

380 

380 

16 

22 

9 

10 

27 

27 

11 

11 

40d. 

781 

781 

26 

40 

9 

9 

33 

34 

13 

13 

40e. 

- 

- 

90 

146 

10 

11 

70 

72 

45 

45 

40f. 

- 

- 

180 

272 

13 

14 

92 

92 

49 

49 

40g. 

- 

- 

206 

319 

23 

25 

123* 

123' 

85 

85 

41a. 

5 

5 

4 

4 

4 

4 

8 

8 

4 

4 

41b. 

6 

6 

4 

5 

4 

5 

18 

18 

4 

4 

41c. 

12 

12 

6 

8 

6 

6 

21 

21 

5 

5 

41d. 

30 

30 

15 

22 

9 

11 

38 

38 

9 

9 

41e. 

222 

222 

29 

38 

17 

20 

47 

47 

14 

14 

41f. 

933 

933 

57 

89 

24 

27 

54 

54 

16 

16 

41g. 

3285 

3285 

84 

144 

29 

30 

62 

62 

21 

21 

32 


Summary  of  Results  :  Nonlinear  Least-Squares  Methods 

(number  of  function  evaluations) 


DeVilliers  and  Glasser  Test  Problems 

DMNG  NPSOL  DMNH  MNA 


42a.° 

4 

24 

53 

56 

60* 

61* 

28 

28 

41* 

41* 

42b.° 

4 

24 

103* 

104* 

140* 

141* 

35 

36 

16 

16 

42c.° 

4 

24 

76 

78 

51* 

52* 

30 

31 

6 

6 

42d.° 

4 

24 

61 

64 

56* 

57* 

30 

30 

6 

6 

43a.° 

5 

16 

49 

51 

44*' 

53*' 

22 

22 

30* 

30* 

43b.° 

5 

16 

58 

60 

37*' 

38" 

26 

27 

17* 

17* 

43c ,° 

5 

16 

41 

44 

44*' 

54*' 

21 

21 

89* 

89* 

43d.° 

5 

16 

57 

60 

112*' 

120*' 

27* 

28* 

41* 

42* 

43e.° 

5 

16 

51 

53 

95* 

97* 

28* 

29* 

142* 

143* 

43f.° 

5 

16 

45 

48 

56* 

59* 

17 

18 

37* 

37* 

Dennis,  Gay,  and  Vu  Test  Problems 


n  m  DMNG  NPSOL  DMNH  MNA 


44a.° 

6 

6 

441 

444 

488 

490 

179 

180 

143 

144 

44b. 0 

6 

6 

31 

34 

57 

59 

9 

10 

46 

46 

44c. 0 

6 

6 

3726 

3731 

- 

- 

194 

195 

914 

915 

44d.° 

6 

6 

_t 

3865 

- 

- 

187 

188 

915 

916 

44e.° 

6 

6 

_t 

2815 

1976 

1978 

219 

220 

475 

476 

45a.° 

8 

8 

284 

288 

474 

476 

63 

64 

186 

186 

45b.° 

8 

8 

36 

40 

82 

84 

15 

16 

38 

38 

45c. 0 

8 

8 

6197 

6200 

- 

- 

321 

322 

1416 

1416 

45d.° 

8 

8 

7929 

7934 

1654 

1656 

328 

329 

1478 

1479 

45e.° 

8 

8 

3341 

3346 

- 

- 

351 

352 

1441 

1441 

33 


Summary  of  Results  :  Unconstrained  Optimization  Methods 

(number  of  function  evaluations) 


DeVilliers  and  Glasser  Test  Problems 
G-N  LMDER  DN2G  LSQFDQ  LSQSDN 


42a.° 

51* 

51* 

18 

19 

29* 

29* 

73*' 

73*' 

60*' 

64*' 

42b.° 

611* 

- 

48* 

49* 

74* 

74* 

94* 

94* 

75’ 

75* 

42c.° 

27* 

27* 

20* 

20* 

32* 

32* 

46* 

46* 

48* 

48* 

42d.° 

24* 

24* 

15* 

16* 

23 

24 

27* 

27* 

27* 

27* 

43a.° 

23* 

23* 

14* 

15* 

31 

32 

33* 

33* 

25* 

33* 

43b.° 

15*' 

20* 

18* 

18* 

20 

20 

45* 

45* 

37* 

45* 

43c.° 

24* 

24* 

11* 

11* 

34' 

41' 

33* 

33* 

25* 

33* 

43d. 0 

18* 

18* 

22* 

23* 

17 

17 

38* 

38* 

30* 

38* 

43e.° 

31* 

31* 

12* 

13* 

28 

29 

27* 

27* 

19* 

27* 

43f.° 

22 

22 

12 

13 

20 

20 

31 

31 

23 

31 

Dennis,  Gay,  and  Vu  Test  Problems 
G-N  LMDER  DN2G  LSQFDQ  LSQSDN 


44a.° 

171 

171 

37 

38 

58 

59 

97 

97 

93 

95 

44b.° 

5 

5 

5 

6 

7 

7 

10 

10 

6 

10 

44c.° 

55 

55 

108 

109 

93 

94 

47 

47 

47 

47 

44d.° 

35 

35 

98 

99 

97 

98 

40 

40 

40 

40 

44e.° 

42 

42 

82 

83 

83 

83 

47 

47 

47 

47 

45a.° 

171 

171 

47 

48 

65 

66 

97 

97 

93 

95 

45b.° 

5 

5 

5 

6 

8 

8 

12 

12 

6 

12 

45c. 0 

41 

41 

164 

165 

129 

130 

47 

47 

47 

47 

45d.° 

35 

35 

144 

145 

168 

168 

42 

42 

42 

42 

45e.° 

42 

42 

130 

131 

173 

173 

49 

49 

43 

43 

34 


5.  Appendix  :  Tables  of  Numerical  Results 


In  this  section,  numerical  results  are  presented  for  a  large  set  of  test  problems,  using 
software  based  on  the  unconstrained  optimization  techniques  and  methods  for  nonlinear 
least  squares  problems  discussed  in  Sections  2  and  3. 


5.1  Sources  and  Presentation 


The  following  is  a 

list  of  the  software  sources 

that  were  used  to  obtain  the  results 

subroutine 

source 

problem  type 

derivatives 

DMNG/SUMSOL 

PORT 

unconstrained 

first 

NPSOL 

SOL  /  NAG 

unconstrained 

first 

DMNH/HUMSOL 

PORT 

unconstrained 

second 

MNA 

NPL  /  NAG 

unconstrained 

second 

G-N 

uses  SOL  /  NAG  LSSOL 

least  squares 

first 

LMDER 

MINPACK 

least  squares 

first 

DN2G/NL2S0L 

PORT 

least  squares 

first 

LSQFDQ 

NPL  /  NAG 

least  squares 

first 

LSQSDN 

NPL  /  NAG 

least  squares 

second 

ACM  -  Association  for  Computing  Machinery 
MINPACK  -  Argonne  National  Laboratory,  U.  S.  A. 

NAG  -  Numerical  Algorithms  Group 

NPL  •  National  Physical  Laboratory,  England 

PORT  -  PORT  Mathematical  Software  Library,  A.  T.  &  T.  Bell  Laboratories,  Inc. 
SOL  -  Systems  Optimization  Laboratory,  Stanford  University 

All  of  the  programs  were  run  in  FORTRAN  using  double  precision  on  the  IBM  3081  and 
IBM  3033  computers  at  Stanford  Linear  Accelerator  Center,  for  which 

relative  machine  precision  eM  =  2.22...  x  10“ 16 ;  y/T^  =  1.49...  x  10-8. 

In  the  tables,  associated  with  the  label  'est.  err.’,  we  include  the  quantity 


where  /*  is  the  value  of  /  at  the  point  of  termination,  and  H/beatllj  is  the  best  available 
estimate  of  the  norm  of  the  solution,  in  order  to  get  some  idea  of  the  error  in  ||/*||2. 
For  those  problems  that  have  nonzero  residuals,  the  value  of  ||/be3t  II2  >s  given  to  six 
figures  of  accuracy,  rounded  down. 

The  following  abbreviations  are  used  in  the  headings  of  the  tables : 

est.  err.  -  error  estimate  (5.1.1) 

conv.  -  termination  conditions 

The  following  abbreviations  are  used  in  the  tables  to  describe  conditions  under  which 
the  algorithm  terminates  abnormally: 

f  lim.  -  function  evaluation  limit  reached 

time  -  time  limit  exceeded 

loop  -  subroutine  appears  to  loop 

For  information  on  the  test  problems,  see  Section  5.9. 
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5.2  Trust-Region  Methods  for  Unconstrained  Optimization 
(PORT/ ACM  DMNH/HUMSL  and  DMNG/SUMSL) 

5.2.1  Software  and  Algorithms 

The  results  were  obtained  using  subroutines  DMNH  and  DHNG,  which  are  double  precision 
versions  of  the  ACM  algorithms  HUMSL  and  SUMSL  available  in  the  PORT  Library  [1984].  A 
subproblem  of  the  form 

min  Qk{p)  =  gjp  +  tj  PTHkp 
subject  to||£>tp||2  <  6k 

is  solved  at  each  iteration  for  the  step  pk  to  the  next  iterate,  where  Dk  is  a  diagonal  scaling 
matrix,  and  Hk  is  the  exact  Hessian  matrix  at  xk  in  DMNH,  and  a  quasi- Newton  approximation 
in  DMNG  (see  Sections  2.3.2  and  2.4). 

5.2.2  Parameters 

Parameters  were  kept  at  their  default  values  with  the  following  exceptions : 

IV(MXFCAL)  -  min  {9999,  lOOOn}  function  evaluation  limit 

IV(MXITER)  -  min  {9999,  lOOOn}  iteration  limit 

V(AFCTOL)  -  TOL2  (varied;  see  tables)  absolute  function  convergence  tolerance 
V(RFCTOL)  -  TOL  (varied;  see  tables)  relative  function  convergence  tolerance 

V(  SCTOL)  -  eM  singular  convergence  tolerance 

V (  XCTOL)  -  TOL  (varied;  see  tables)  x  convergence  tolerance 

V(  XFTOL)  -  eM  false  convergence  tolerance 

V(  LMAXO)  -  usually  1.0  (default)  t  initial  trust-region  diameter 

f  In  some  cases  the  default  V(LMAXO)  =  1.0  for  the  initial  diameter  of  the  trust-region  was  too 
large  and  overflow  occurred  during  function  evaluation.  These  cases  are  indicated  in  the  table 
by  giving  the  lower  value  of  V(LHAXO)  that  was  subsequently  used  to  obtain  the  results  in  the 
column  labeled  “init.  diam.”. 

See  Dennis,  Ga^,  and  Welsch  [1981a,  1981b],  Gay  [1983],  and  PORT  [1984]  for  details 
concerning  the  parameters. 
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5.2.3  Convergence  Criteria 


The  following  quantities  will  be  used  in  describing  the  convergence  criteria : 


objective  function 

T*  (=UT/i) 

objective  gradient 

9k  =  VTk  (=  J?fk) 

current  step 

Pi,  the  minimi zer  of  the  subproblem 

Newton  step 

f  H^g  if  Hk  is  positive  definite ; 

\  undefined  otherwise. 

Newton  reduction 

_  f  -Qt(pn)  if  Hk  is  positive  definite ; 
N  to  otherwise. 

predicted  reduction 

Pp  =  ~Qk(Pk) 

actual  reduction 

Pa  =  xk  +Pl) 

scaled  distance 

lYr  ,,  m  maxi<<<"  {\{D(x  -  y))i\} 

1  ,V'  }  ~  ffiaxKf<n  (\{Dxi\+  |(Hy),|}’  T 

t  Here  v,-  denotes  the  ith  component  of  the  vector  v.  There  is  a  provision  for  the  user  to  replace 
the  function  v\  we  used  the  default  in  all  of  the  tests. 

The  convergence  criteria  used  in  DMNH  and  DKNG  are  as  follows : 

•  Absolute  function  convergence  occurs  at  x*  if 

\Tk\  <  V(AFCTOL).  (5.2.1) 

•  Relative  function  convergence  is  intended  to  approximate  the  condition 

Tk  -  T{x')  <  V(RFCTOL)  \Fk\- 


The  test  actually  used  is 


pN  <V(RFCT0L)  |*fc|. 


•  x  convergence  is  intended  to  approximate  the  condition 


(5.2.2) 


*(**,*•,£>*)  <V(XCT0L), 


The  test  actually  used  is 


pt  =  Pn  and  *(*i,*i+pi,Di)<V(XCTOL).  (5.2.3) 

•  Singular  convergence  is  intended  to  approximate  the  condition 

Tk  -  min  (:F(y)  |  ||D*(y  -  *fc)ll  <  V(LMAXS)}  <  V(SCTOL)  \Tk\, 
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where  Dk  is  the  diagonal  scaling  matrix  at  the  ifcth  iterate.  The  test  for  singular  convergence 
is  made  only  when  none  of  the  convergence  criteria  listed  above  holds.  It  is  meant  to  indicate 
relative  function  convegence  when  the  Hessian  in  the  subproblem  is  singular. 

The  actual  test  is 

Tk  -  min  {<2fc(y)  |  ||Dk(y  -  c*)ll  <  V(LMAXS)}  <  V(SCTOL)  \Tk\ .  (5.2.4) 

Under  certain  conditions,  the  test  (5.2.4)  is  repeated  for  a  step  of  length  V(LMAXS). 

•  False  convergence  is  returned  if  none  of  the  other  criteria  are  satisfied  and  a  trial  step  no 
larger  than  V(XFCTOL)  is  rejected.  This  usually  indicates  either  an  error  in  computing  the 
objective  gradient,  a  discontinuity  (in  T  or  g)  near  the  current  iterate,  or  that  one  or  more  of 
the  convergence  tolerances  (V(RFCTOL),  V(XCTOL),  and  V(AFCTOL))  are  too  small  relative 
to  the  accuracy  to  which  the  objective  is  computed. 

The  test  actually  used  is 

Tk  —  F(xk  +  pk)  <  V(TUNERl)pP  and  v(xk,xk  +pk>Dk)  <  V(XFTOL),  (5.2.5) 

where  the  parameter  V(TUNERl)  is  adjustable,  although  in  these  tests  the  default  value  0.1 
is  used  throughout. 

Except  for  (5.2.1),  tests  for  convergence  are  performed  only  when 

Pa  <  2 pP.  (5.2.6) 

See  Dennis.  Gay,  and  Welsch  (1981a,  1981b],  Gay  [1983],  and  PORT  [1984]  and  PORT 
[1984]  for  more  discussion  of  the  convergence  criteria. 

The  following  abbreviations  are  used  in  the  tables  to  describe  the  conditions  under  which  the 
algorithm  terminates : 


ABS  F  - 

(5.2.1) 

REL.  F  - 

(5.2.2)  and  (5.2.6) 

X 

(5.2.3)  and  (5.2.6) 

X,  F 

(5.2.2)  and  (5.2.3)  and  (5.2.6) 

SING  - 

(5.2.4)  and  (5.2.6) 

FALSE  - 

(5.2.5)  and  (5.2.6) 

F  LIM  - 

function  evaluation  limit  reached 

TIME  - 

time  limit  exceeded 

LOOP  - 

subroutine  appears  to  loop 

39 


The  total  number  of  Jacobian  evaluations  is  either  equal  to  the  total  number  of  iterations 
of  the  method,  or  it  is  one  more  than  the  number  of  iterations.  The  number  in  the  column 
labeled  “iters.  /  J  evals.”  is  followed  by  a  “+"  if  an  extra  Jacobian  evaluation  was  used  in 
the  computation. 
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5.3  Quasi-Newton  (BFGS)  Linesearch  Method 

(SOL/NAG  NPSOL) 

5.3.1  Software  and  Algorithm 

The  results  were  obtained  using  subroutine  NPSOL  from  the  Systems  Optimization  Lab¬ 
oratory  (Systems  Optimization  Laboratory,  Department  of  Operations  Research,  Stanford 
University),  Stanford  University,  also  available  in  the  NAG  Library.  In  NPSOL  a  search  direc¬ 
tion  is  determined  at  each  iteration  from  a  subproblem  of  the  form 

x  3k  P  +  ^  PTHkP, 

where  the  Hessian  matrix  Hk  is  calculated  using  the  BFGS  method  initialized  with  I  (see  Sec¬ 
tion  2.4.2).  This  is  followed  by  a  linesearch  that  uses  both  function  and  gradient  information 
to  obtain  a  steplength  along  the  search  direction  [Gill  et  al.  (1979)]. 

5.3.2  Parameters 

Parameters  were  kept  at  their  default  values  with  the  following  exceptions  f  : 

Infinite  Bound  Size  -  1020 

Infinite  Step  Size  -  1020 

Iteration  Limit  -  lOOOOf 

Optimality  Tolerance  -  varied;  see  tables 
Step  Limit  -  usually  2.0  (default)  1 

f  In  those  cases  (Problems  44c.,  d.  (n  =  m  =  6)  and  45c.,  e.  (n  =  m  =  8)  in  which  the 
iteration  limit  was  actually  reached,  the  results  listed  in  the  tables  are  taken  from  first  iteration 
in  which  the  number  of  function  evaulations  reaches  or  exceeds  lOOOn. 

1  In  some  cases  the  default  Step  Limit  =  2.0  was  too  large  and  overflow  occured  during  function 
evaluation  in  the  linesearch.  These  cases  are  indicated  in  the  tables  by  giving  the  lower  value  of 
Step  Limit  that  was  subsequently  used  to  obtain  the  results  in  the  column  labeled  “Step  Lim” . 

See  Gill  et  al.  [1986]  for  details  concerning  the  parameters. 

5.3.3  Convergence  Criteria 
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The  following  quantities  will  be  used  in  describing  the  convergence  criteria  : 

objective  function  (=  1  /*) 

objective  gradient  gk  =  VJ7*  (=  Jj ft) 

optimality  tolerance  fop* 

The  sequence  of  iterates  generated  by  NPSOL  is  judged  to  have  converged  if  the  following 
two  conditions  hold : 


<**IWI2  <  +  ||**||2) 

(5.3.1) 

and 

IM2  <  y/*opt(\  +  max  {(1  +  |^t|)||fft||2}) 

(5.3.2) 

M  <  (1  +  max  {(1  +  |**|),||*||2}). 

(5.3.3) 

Condition  (5.3.1)  is  meant  to  ensure  that  the  sequence  {z*}  has  converged,  while  conditions 
(5.3.2)  and  (5.3.3)  are  intended  to  test  whether  the  requirement  that  the  gradient  vanish  is 
approximately  satisfied  at  xk.  Condition  (5.3.3)  allows  NPSOL  to  accept  a  point  as  a  local 
mimimum  if  a  more  restrictive  test  on  the  necessary  condition  than  (5.3.2)  is  satisfied,  but 
condition  (5.3.1)  does  not  hold.  For  a  detailed  discussion  of  convergence  criteria  similar  to 
these,  see  Section  8.2  of  Gill,  Murray,  and  Wright  [1981]. 

The  following  abbreviations  are  used  in  the  tables  to  describe  the  conditions  under  which  the 
algorithm  terminates  :  f 

opt.  -  optimal  point  found 
*  -  current  point  cannot  be  improved 

**  -  optimal  solution  found,  but  requested  accuracy  could  not  be  achieved 

f  lim.  -  function  evaluation  limit  reached 

t  A  corresponds  to  the  situation  in  which  the  algorithm  terminates  due  to  failure  in  the 
linesearch  to  find  an  acceptable  step  at  the  current  iteration.  A  occurs  when  condition 
(5.3.1)  is  satisfied  but  not  condition  (5.3.2) ;  that  is,  conditions  for  optimality  are  met  at  the 
current  point  but  the  iterates  have  not  yet  converged. 
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Numerical  Results  for  HPSOL 
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5.4  Second-Derivative  (Modifled-Newton)  Linesearch  Method 

(NPL/NAG  MNA) 

5.4.1  Software  and  Algorithm 

The  results  were  obtained  using  subroutine  MNA  from  the  National  Physical  Laboratory, 
available  at  Stanford  Linear  Accelerator  Center.  The  algorithm  implements  a  modified  New¬ 
ton  method  in  which  the  search  direction  at  each  iteration  is  the  solution  to  a  subproblem 
of  the  form 

min  ffitP  +  | PTRkP, 

and  the  exact  Hessian  matrix  is  replaced  by  modified  Cholesky  factors  if  it  is  either  indefinite 
or  computationally  singular  (see  Gill  and  Murray  [1974a]  and  Section  2.4.1).  A  step  length 
along  the  search  direction  is  then  computed  by  a  linesearch  method  [Gill  and  Murray  (1974b)] 
that  uses  both  function  and  gradient  information  to  obtain  sufficient  decrease  in  the  objective 
function.  MNA  requires  exact  second  derivatives,  and  is  similar  to  subroutine  E04LBF  from 
the  NAG  Library  [1984],  the  principal  difference  being  that  the  latter  allows  specification  of 
fixed  upper  and  lower  bounds  on  the  variables. 

5.4.2  Parameters 

Parameters  were  kept  at  their  default  values  with  the  following  exceptions  : 

MAXCAL  -  min  {9999,  lOOOn}  function  evaluation  limit 

XTOL  -  varied;  see  tables  accuracy  in  x 

ETA  -  0.9  linesearch  accuracy 

STEPMX  -  usually  10®  (default)  f  maximum  step  for  linesearch 

f  In  some  cases  the  default  STEPMX  =  106  was  too  large  and  overflow  occurred  during  function 
evaluation  in  the  linesearch.  These  cases  are  indicated  in  the  table  by  giving  the  lower  value  of 
STEPMX  that  was  subsequently  used  to  obtain  the  results  in  the  column  labeled  “max.  step” . 

See  NAG  [1984]  for  details  concerning  the  parameters. 
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5.4.3  Convergence  Criteria 

The  following  quantities  will  be  used  in  describing  the  comvergence  criteria : 

objective  function  (=  ^  fj ft) 

objective  gradient  :  gt  =  VTk  (=  ft) 
search  direction  Pk,  the  minimizer  of  the  subproblem 

steplength  a determined  by  the  linesearch 

An  iterate  is  determined  to  be  optimal  by  MNA  if  the  following  four  conditions  hold  : 

a*!|p*||2(XT0L  +  y/eZKl  +  ||zt||2)  (5.4.1) 

and 

Fk- 1  -  Fk  <  (XTOL2  +  cM)(l  +  |JTt|)  (5.4.2) 

and 

||<,fc||2<(XT0L  +  ci/3)(l  +  in|)  (5.4.3) 

and 


or  if 


V2^t  is  positive  definite, 

(5.4.4) 

lb*!!’  <  Q-Oly/eZ- 

(5.4.5) 

A  necessary  condition  for  optimality  is  that  the  gradient  vanish,  and  conditions  (5.4.3)  and 
(5.4.5)  are  intended  to  test  whether  this  requirement  is  approximately  satisfied  at  z*.  Con¬ 
ditions  (5.4.1)  and  (5.4.2)  are  meant  to  ensure  that  the  sequence  {z*}  has  converged,  while 
condition  (5.4.4),  together  with  condition  (5.4.3),  implies  that  sufficient  conditions  for  a 
strict  local  minimum  appear  to  hold  at  z*.  Condition  (5.4.5)  allows  MNA  to  accept  a  point  as 
a  local  mimimum  if  a  more  restrictive  test  than  (5.4.1)  on  the  necessary  condition  is  met,  but 
one  or  more  of  the  other  conditions  for  convergence  do  not  hold.  For  a  detailed  discussion 
of  convergence  criteria  similar  to  these,  see  Section  8.2  of  Gill,  Murray,  and  Wright  [1981]. 


The  following  abbreviations  are  used  in  the  tables  to  describe  the  conditions  under  which  the 
algorithm  terminates : 

opt.  -  optimal  point  found 

*  -  current  point  cannot  be  improved  f 

F  um.  -  function  evaluation  limit  reached 
time  -  time  limit  exceeded 

t  A  V  corresponds  to  the  situation  in  which  the  algorithm  terminates  due  to  failure  in  the 
linesearch  to  find  an  acceptable  step  at  the  current  iteration. 
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5.5  Gauss-Newton  Methods 


5.5.1  Software  and  Algorithm 

The  software  package  LSSOL  [Gill  et  al.  (1986a)]  is  used  to  solve  the  linear  least-squares 
subproblem  (3.1.3).  The  linesearch  procedure  used  for  the  numerical  examples  in  this  section, 
requires  both  function  and  gradient  information.  It  is  taken  from  the  nonlinear  programming 
code  NPSOL  [Gill  et  al.  (1979);  (1986b)]. 


5.5.2  Parameters 


Parameters  in  LSSOL  were  kept  at  their  default  values  with  the  following  exceptions : 

Rank  Tolerance  -  varied,  see  tables 
Infinite  Bound  Size  -  1020 

See  Gill  et  al.  [1986a]  for  details  concerning  the  parameters. 

In  addition,  the  following  parameters  are  chosen  for  the  linesearch  : 

t]  -  0.5 

“max  -  min  {(100(1 +  ||z||2)+ l)/||p||2,1020}  f 
f  In  some  cases  the  default  value  amax  was  too  large  and  overflow  occurred  during  function 
evaluation  in  the  linesearch.  These  cases  are  indicated  in  the  tables  by  giving  the  value  y  <  100 
such  that  ama>:  ==  min{(7(l  +  ||z||2)  +  l)/|jp||2, 1020)}  that  was  subsequently  used  to  obtain 
the  results  in  the  column  labeled  “step  fac.”. 

See,  e.  g.,  Gill,  Murray,  and  Wright  [1981]  for  a  discussion  of  the  linesearch  parameters. 

5.5.3  Convergence  Criteria 

Convergence  is  judged  to  have  occurred  at  the  ifcth  iterate  if  either 

IIAII,  <  <°«9  (5.5.1) 

or 


MI2<eL/3(i  +  IIAII2). 


(5.5.2) 
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The  algorithm  is  also  terminated  if  there  is  a  negligible  change  in  x, 

«fcl|P*lla<4i*(l  +  IM,).  (5-5.3) 

where  a *  is  the  step  length  determined  by  the  linesearch. 

5.5.4  Table  Information 

Under  the  label  ‘conv.’,  the  following  notation  is  used  to  describe  conditions  under  which 
the  algorithm  terminates  : 

ABS.  F  -  (5.5.1) 

g  -  (5.5.2) 

x  -  (5.5.3) 

f  lim.  -  function  evaluation  limit  reached 

A  superscript  0  following  a  problem  number  indicates  a  zero-residual  problem. 

A  superscript  L  following  a  problem  number  denotes  a  linear  least-squares  problem. 
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10  13 

10-3° 

ABS.  P,  Q 

26a.° 

10 

10 

1.49x10  -8 

8 

.306 

10"11 

10"u 

10-22 

O 

2.23x10  -16 

— 

8 

.306 

10"n 

10- 11 

10-22 

O 

26b.° 

20 

20 

1.49x10  -8 

mm 

pa 

10-14 

10-28 

O 

2.23x10  -16 

mm 

mm 

10-28 

Q 

27a.° 

10 

10 

1.49x10  -8 

mm 

mm 

3.18 

10-15 

10-14 

10"29 

ABS.  P,  a 

2  23x10  -16 

H 

KM 

3.18 

10-15 

10-14 

10” 29 

ABS.  P.  a 

27b.° 

20 

20 

1.49x10  -8 

10.0 

31 

KM 

Q 

2.23x10  -16 

10.0 

31 

Hi 

O 

28a.° 

10 

10 

1.49x10  -8 

KM 

.412 

10-15 

10-16 

10-si 

ABS.  P,  O 

2.23x10  -16 

.412 

10  1 5 

10-16 

10-51 

ABS.  P,  O 

28b.° 

20 

20 

1  49x10  -8 

4 

3 

.571 

10-16 

10-16 

10"32 

ABS.  P,  O 

2  23x10  -16 

4 

3 

.571 

10-16 

10-16 

10-32 

ABS.  P,  a 
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n 

m 

rank 

tol. 

step 

fac. 

f,J 

evals. 

iters. 

11*% 

uni. 

11*% 

est. 

err. 

conv. 

29a.° 

10 

10 

1.49x10  -8 

mm 

.412 

10"14 

10-14 

10-29 

ABS.  P,  a 

2  23x10  -IS 

*mm 

.412 

10~14 

10“ 14 

10-29 

ABS.  F,  a 

29b.° 

20 

20 

1.49x10  -8 

4 

3 

.571 

a 

2.23x10  -16 

4 

3 

.571 

a 

30a.° 

10 

10 

1.49x10  -8 

6 

5 

2.05 

10-16 

10-15 

10-3! 

ABS.  P,  0 

2.23x10  -16 

6 

5 

2.05 

IQ-16 

10-15 

10-3! 

ABS.  P.  0 

30b.° 

20 

20 

1.49x10  -8 

6 

5 

3.04 

10-16 

10-15 

IO"31 

ABS.  P,  0 

2.23x10  -16 

6 

5 

3.04 

10-16 

IQ-15 

10-3! 

ABS.  P,  0 

31a.° 

10 

10 

1.49x10  -8 

7 

6 

1.80 

IQ-15 

10-is 

10-31 

ABS.  P,  ◦ 

2.23x10  -16 

7 

6 

1.80 

IO"15 

10-15 

IO"31 

ABS.  P,  Q 

31b.° 

20 

20 

1.49x10  -8 

7 

6 

2.66 

10-!5 

10-15 

10-si 

ABS.  P,  a 

2.23X10  -16 

7 

6 

2.66 

IQ-15 

IQ-15 

IO-31 

ABS.  P,  O 

32/ 

10 

20 

1.49x10  -8 

2 

1 

3.16 

10° 

IO"14 

0.00 

a 

2.23X10  -16 

2 

1 

3.16 

10° 

IO"14 

0.00 

0 

33/ 

10 

20 

1.49x10  -8 

3 

2 

Hi 

10~10 

a,  x 

2.23x10  -16 

8 

8 

ES9 

1TM 

gTP  >  0 

34/ 

10 

20 

1.49x10  -8 

3 

2 

4.90 

10° 

10~u 

IO-6 

a,  x 

2.23x10  -16 

3 

2 

4.90 

10° 

io-u 

IO"6 

o,  X 

35&i 

8 

8 

1.49X10  -8 

3053 

386 

1.61 

10"1 

io- 1 

IO"3 

X 

2  23x10  -16 

(8003) 

(1012) 

1.60 

10-1 

10° 

IO"2 

P  LIM. 

^5b.° 

9 

9 

1.49x10  -8 

148 

29 

1.73 

10~14 

10-14 

IO-29 

ABS.  P,  a 

2  23x10  -16 

249 

38 

1.70 

10- 1 

10° 

IO-2 

X 

35c. 

10 

10 

1.49x10  -8 

(10006) 

(1353) 

1.79 

10-1 

10° 

IO"2 

P  LIM. 

2  23x10  -16 

232 

31 

1.79 

10-1 

10° 

IO"2 

X 

36a.° 

4 

4 

1.49x10  -8 

490 

18.8 

10-6 

10-12 

IO-12 

O 

2.23x10  -16 

21 

50.0 

10-15 

10~14 

IO"31 

ABS.  P,  O 

36b.° 

9 

9 

1.49x10  -8 

683 

128 

18.8 

10~6 

IO"12 

10~12 

a 

2  23x10  -16 

36 

21 

50.0 

10-15 

10" 14 

IO-31 

ABS.  P,  O 

36c. 0 

9 

9 

1.49x10  -8 

20 

19 

1.73 

10'11 

IO"11 

10-22 

O 

2.23x10  -16 

20 

19 

1.73 

nr11 

IO"11 

10-22 

O 

36d.° 

9 

9 

1.49x10  -8 

74 

29 

19.1 

IO-6 

IO-11 

10-12 

O 

2.23x10  -16 

(9001) 

(1190) 

345. 

10"6 

10-5 

10-13 

P  LIM. 

37. 

2 

16 

1  49x10  -8 

39 

38 

8.85 

101 

IO-8 

10~6 

X 

2.23x10  -16 

39 

38 

8.85 

101 

10“8 

10“6 

X 

38. 

3 

16 

1.49X10  -8 

58 

56 

26.1 

101 

10- 10 

IO"6 

0 

2  23x10  -16 

58 

56 

26.1 

101 

IO-10 

IO-6 

0 
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n 

m 

rank 

step 

/,  J 

iters. 

11**112 

lint 

Ib'lla 

est. 

conv. 

tol. 

fac. 

evals. 

err. 

39a. 

2 

3 

1.49X10  -8 

8 

7 

10"® 

10"1 

10~n 

io-7 

a 

2.23X10  -16 

8 

7 

10"® 

10"1 

10~n 

io-7 

a 

39b. 

2 

3 

1.49X10  -8 

IO"7 

10"1 

10~n 

io-7 

a 

2.23X10  -16 

IO"7 

10-1 

10-11 

IO-7 

a 

39c. 

2 

3 

1.49X10  -8 

23 

14 

10~7 

10"1 

10-12 

IO"7 

a 

2.23X10  -16 

23 

14 

10- 7 

10"1 

10-12 

io-7 

a 

39d. 

2 

3 

1.49X10  -8 

681 

333 

10"7 

10-1 

10-10 

io-7 

a 

2.23x10  -16 

681 

333 

10"7 

10"1 

10-1° 

IO"7 

o 

39e. 

2 

3 

1.49X10  -8 

(2001) 

(925) 

10"7 

10_1 

10~9 

io-7 

P  L1M. 

2.23X10  -16 

(2001) 

(925) 

10"7 

IO"1 

10~9 

10-7 

P  LIM. 

39f. 

2 

3 

1.49X10  -8 

(2000) 

(873) 

10"9 

IO-1 

10~8 

10-7 

P  LIM. 

2  23X10  -16 

(2000) 

(873) 

10"9 

10"1 

00 

1 

o 

10-7 

P  LIM. 

39g. 

2 

3 

1.49X10  -8 

(2001) 

(607) 

10~8 

10"1 

10-s 

IO-7 

P  LIM. 

2.23x10  -16 

(2001) 

(607) 

10"8 

10"1 

10~® 

10-7 

P  LIM. 

40a. 

3 

4 

1.49x10  -8 

13 

12 

10~6 

10° 

10“n 

10-7 

a 

2  23x10  -16 

13 

12 

10“6 

10° 

10-11 

IO-7 

O 

40b. 

3 

4 

1.49X10  -8 

16 

10 

10“6 

10° 

10-12 

IO-7 

O 

2.23x10  -16 

16 

10 

10~® 

10° 

10-12 

IO"7 

a 

40c. 

3 

4 

1  49x10  -8 

188 

10"7 

10° 

10-10 

IO"7 

a 

2.23x10  -16 

188 

10-7 

10° 

10-10 

IO-7 

o 

40d. 

3 

4 

1.49X10  -8 

781 

ESI 

10-7 

10° 

10-10 

IO-7 

o 

2.23x10  -16 

781 

ESI 

10-7 

10° 

10~1C 

IO-7 

a 

40e. 

3 

4 

1.49x10  -8 

(3000) 

(690) 

o 

1 

to 

o 

o 

10"1 

CO 

1 

o 

P  LIM. 

2.23x10  -16 

(3000) 

(690) 

o 

\ 

to 

10° 

10- 1 

10~3 

P  LIM. 

40f. 

3 

4 

1.49x10  -8 

(3002) 

(630) 

.101 

o 

o 

101 

IO"1 

P  LIM. 

2.23x10  -16 

(3002) 

(630) 

.101 

o 

o 

101 

IO"1 

P  LIM. 

40g. 

3 

4 

1  49X10  -8 

(3001) 

(607) 

.110 

10l 

103 

101 

P  LIM. 

2.23x10  -16 

(3001) 

(607) 

.110 

101 

103 

IO1 

P  LIM. 

41a. 

5 

10 

1.49X10  -8 

5 

4 

10-6 

10° 

10-18 

IO-7 

a 

2.23X10  -16 

5 

4 

10"6 

10° 

10~13 

IO-7 

a 

41b. 

5 

10 

1  49X10  -8 

10-6 

10° 

10-1° 

IO-7 

a 

2.23x10  -16 

10"® 

10° 

10-10 

IO"7 

o 

41c. 

5 

10 

1.49x10  -8 

12 

11 

10-6 

10° 

10-11 

IO"7 

a 

2.23x10  -16 

12 

11 

10"® 

10° 

10~u 

IO-7 

a 

41d.  5  10  1 49xio  -8  30  20  10-®  10°  IO-10  10-7 

2.23X10  -16  30  20  10"®  10°  10~10  10-7 


41e.  5  10  1 49xio  -s  222  110  10"7  10°  10~10  10-7 

_ 2.23X10  -16 _ 222  110  IQ-7  10°  1Q~10  10~7 

41f.  5  1  0  1.49X10  -8  933  414  10-7  10°  IO-10  IO-7 

_ 2.23x10  -16 _ 933  414  IQ"7  10°  1Q~10  1Q~7 

41g.  5  10  1 49xio  -8  3285  1065  IO-8  10°  IO-10  IO"7 

2  23X10  -16  3285  1065  IO-8  10°  IO-10  IO-7 
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n 

m 

rank 

tol. 

step 

fac. 

/,  J 

eva la. 

iters. 

II**  1 1 2 

Il/*ll2 

n 

est. 

err. 

conv. 

42a.° 

4 

24 

1.49x10  -8 

0.1 

mm 

— 

2.23x10  -16 

0.1 

mm. 

■H 

42b. 0 

4 

24 

1.49x10  -8 

0.1 

611 

316 

61.9 

IQ-13 

o 

o 

*— 4 

io- 25 

X 

2.23x10  -16 

0.1 

(4002) 

(1577) 

105 

102 

10u 

105 

P  L1M. 

42c. 0 

4 

24 

1.49x10  -8 

0.1 

mm 

mm 

60.3 

10-13 

10"n 

10-26 

Q 

2.23x10  -16 

0.1 

mm 

60.3 

io- 13 

IO"11 

IO-26 

o 

42d.° 

4 

24 

1.49x10  -8 

0.1 

24 

23 

60.3 

io- 13 

10"n 

10-26 

o,  X 

2.23x10  -16 

0.1 

24 

23 

60.3 

io-13 

10-11 

10-26 

o,  X 

43a.° 

5 

16 

1.49x10  -8 

1.0 

23 

mm 

54.0 

10"14 

10-12 

10-28 

o 

2.23x10  -16 

1.0 

23 

mm 

54.0 

IO"14 

10- 12 

IO"28 

o 

43b. 0 

5 

16 

1.49x10  -8 

0.1 

15 

13 

53.6 

10_1 

10-7 

IO-2 

X 

2.23x10  -16 

0.1 

20 

14 

54.0 

io-14 

IO"12 

10-28 

o 

43c. 0 

5 

16 

1.49X10  -8 

1.0 

24 

14 

54.0 

10" 14 

10-11 

IO-27 

o 

2.23x10  -16 

1.0 

24 

14 

54.0 

10"14 

IO-11 

IO-27 

a 

43d. 0 

5 

16 

1.49x10  -8 

1.0 

18 

10 

54.0 

IO-14 

IO"11 

IO-27 

a 

2  23x10  -16 

1.0 

18 

10 

54.0 

10"14 

10-11 

IO-27 

a 

43e.° 

5 

16 

1.49X10  -8 

1.0 

31 

17 

54.0 

10"14 

10-12 

10~28 

a 

2  23x10  -16 

1.0 

31 

17 

54.0 

10-H 

10-12 

10-28 

o 

43f.° 

5 

16 

1  49X10  -8 

22 

15 

54.0 

10- 14 

10-12 

10- 28 

a 

2.23x10  -16 

22 

15 

54.0 

10"14 

10-12 

10-28 

a 

44a. 0 

6 

6 

1  49x10  -8 

171 

45 

4.06 

10-13 

10"n 

10- 25 

a 

2.23x10  -16 

171 

45 

4.06 

10-13 

10" 11 

10-25 

o 

44b. 0 

6 

6 

1.49x10  -8 

5 

4 

3.52 

10-1* 

10-13 

IO"29 

ABS-  P.  a 

2.23x10  -16 

5 

4 

3.52 

10-15 

10-13 

10-29 

ABS.  F,  o 

44c. 0 

6 

6 

1.49x10  -8 

55 

17 

20.6 

10-M 

10-10 

IO-28 

X 

2  23x10  -16 

55 

17 

20.6 

10-1“ 

10-1° 

IO-28 

X 

44d.° 

6 

6 

1.49x10  -8 

mm 

15.3 

10"14 

10-11 

IO-29 

ABS.  P,  a 

2.23x10  -16 

■9 

15.3 

10-14 

10"  i 1 

IO-29 

ABS.  F,  O 

44e.° 

6 

6 

1.49x10  -8 

42 

mm 

9.27 

10-1“ 

10- 12 

10- 28 

ABS.  P,  O,  X 

2.23x10  -16 

42 

■9 

9.27 

10-14 

10-12 

IO-28 

ABS.  P,  a,  X 

45a. 0 

8 

8 

1.49x10  -8 

171 

45 

4.06 

10-13 

10-H 

10-25 

o 

2.23x10  -16 

171 

45 

4.06 

10-13 

10-“ 

IO”25 

o 

45b. 0 

8 

8 

1.49x10  -8 

5 

4 

3.56 

IO"15 

10- 13 

IO'29 

ABS.  P,  O 

2.23x10  -16 

5 

4 

3.56 

IO-15 

10-13 

10-29 

ABS.  F,  Q 

45c. 0 

8 

8 

1.49x10  -8 

41 

17 

20.6 

IO"14 

IO-11 

IO-29 

ABS.  F,  X 

2.23x10  -16 

41 

17 

20.6 

IO"14 

10“n 

10-29 

ABS.  F,  X 

45d.° 

8 

8 

1  49x10  -8 

mm 

15.3 

10-16 

IO-13 

IO"31 

ABS.  P,  a 

2.23x10  -16 

■9 

15.3 

IO-16 

10-13 

10-51 

ABS.  Pt  O 

45e.° 

8 

8 

1.49x10  -8 

mm 

9.31 

10-15 

10-12 

IO"29 

ABS.  P,  O,  X 

2  23x10  -16 

■a 

9.31 

10-15 

10-12 

IO-29 

ABS.  P,  O,  X 
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5.6  Levenberg-Marquardt  Method 

(MINPACK  LMDER) 

5.6.1  Software  and  Algorithm 

The  results  were  obtained  using  the  MINPACK  subroutine  LMDER,  which  implements  a 
Levenberg-Marquardt  method  using  exact  derivative  information.  A  subproblem  of  the  form 

min  Qk(p)  =  2gjp  + pTJk  Jkp 
pes* 

subject  to||L>ip|[2  <  6t 

is  solved  at  each  iteration  for  the  step  pk  to  the  next  iterate,  where  Dk  is  a  diagonal  scaling 
matrix. 

5.6.2  Parameters 

The  results  were  obtained  using  the  MINPACK  subroutine  LMDER,  with  the  following 
input  parameters  : 


XTOL  - 

varied,  see  tables 

accuracy  in  x 

FTOL  - 

varied,  see  tables 

accuracy  in  sum  of  squares 

GTOL  - 

0.00 

gradient  norm  tolerance 

MAXFEV  - 

min  {9999, 1000  *  n) 

function  evaluation  limit 

MODE  - 

1 

specifies  internal  scaling 

FACTOR  - 

100.  (default) 

initial  step  magnification 

f  In  some  cases  the  default  FACTOR  =  100.0  was  too  large  and  overflow  occurred  during  function 
evaluation.  These  cases  are  indicated  in  the  table  by  giving  the  lower  value  of  FACTOR  that  was 
subsequently  used  to  obtain  the  results. 

For  details  about  these  parameters,  More,  Garbow,  and  Hillstrom  [1980], 
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5.6.3  Convergence  Criteria 


The  following  quantities  will  be  used  in  describing  the  convergence  criteria : 


residual  vector 
ith  residual  gradient 
Jacobian  matrix 
objective  function 
objective  gradient 
current  step 

predicted  reduction 


actual  reduction 


/(**) 

A**) 

T{xk)  =  f(xk)Tf(xk) 

gk  =  VjT(*t)  =  2  J(xk)Tf(xk) 

pk ,  the  minimizer  of  the  subproblem 

1IAIU  -  HA  +  JtPkWi  -Qtjpk) 

Pp~  ll/*ll,  '  ll/fcllz 

ll/(«fc)lla  -  ll/(«*  +  p*)M2  Fjxkj  -  H*k  ±  Pk) 

PA  "  ll/kllz  ”  IIMIz 


Criteria  for  termination  of  LMDER  at  xk  are  as  follows  follows : 

•  T  convergence.  Both  actual  and  predicted  reductions  in  the  sum  of  squares  are  at  most 


FTOL. 


| pA\  <  FTOL  and  pP  <  FTOL  and  pA  <  2 pP 


(5.6.1) 


This  attempts  to  guarantee  that 


IIAII2<(i  +  ftol)U/(z*)||2. 


x  convergence.  Relative  error  between  two  consecutive  iterates  is  at  most  XTOL. 


«t+1  <  XTOL||*t  +pfc| 


(5.6.2) 


This  attempts  to  guarantee  that 


II Dk(xk  -  x*)||2  <  XT0L||£»*(x*)|(2. 

•  The  cosine  of  the  angle  between  fk  and  any  column  of  Jk  is  at  most  GTOL  in  absolute 


|  V<f>i(xk)Tfk  . 

max  I  — — ...  rr~  <  GTOL 

I<»<m  ||V^,(zjb)||2||/t|l2 

This  approximates  the  necessary  condition  g(xk)  =  0. 


(5.6.3) 
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•  FTOL  is  too  small.  No  further  reduction  in  the  sum  of  squares  is  possible. 

\Pa\  <  and  pP  <  eu  and  pA  <  2 pP  (5.6.4) 

•  XTOL  is  too  small.  No  further  improvement  in  the  approximate  solution  xt  is  possible. 

<5*+i  <  f*flN  +Pfc||2  (5.6.5) 


The  convergence  criteria  are  described  in  more  detail  in  More,  Garbow,  and  Hillstrom  [1980]. 
The  following  abbreviations  are  used  in  the  tables  to  describe  the  conditions  under  which  the 
algorithm  terminates: 

-  (5.6.1)  and  (5.6.7) 
x  -  (5.6.2)  and  (5.6.7) 

x.  f  -  (5.6.1)  and  (5.6.2)  and  (5.6.7) 
g  -  (5.6.6)  and  (5.6.7) 

F  lim.  -  function  evaluation  limit  reached 
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Numerical  Results  for  LMDER 
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m 
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max. 

step 

/ 

evals. 

iters./ 

J  evals. 

ll**lla 

lirila 

ll^lla 

est. 

err. 

conv. 

l.° 

2 

2 

IO-8 

22 

16 

1.41 

10_1® 

IO"15 

10~32 

X 

IO-12 

22 

16 

1.41 

10-1® 

10-15 

10-32 

X 

2.° 

2 

2 

10~8 

14 

8 

11.4 

101 

10“3 

101 

REL.  F 
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5.7  Adaptive  Special  Quasi-Newton  Method 

(PORT/ ACM  DN2G/NL2S0L) 


5.7.1  Software  and  Algorithm 

The  results  were  obtained  using  subroutine  DN2G,  a  double  precision  version  of  the 
ACM  algorithm  NL2S0L  available  in  the  PORT  Library  [1984].  A  subproblem  of  the 
form 

min  Qk(p )  =  gjp  +  \pT(JkJk  +  Bk)p 

subject  to||£)fcpj|2  <  6k 

is  solved  at  each  iteration  for  the  step  pk  to  the  next  iterate,  where  Dk  is  a  diagonal 
scaling  matrix.  The  method  is  adaptive,  so  that  Bk  is  sometimes  null  and  sometimes 
a  scaled  quasi-Newton  approximation  to  the  part  of  the  Hessian  involving  the  second 
derivatives  of  /. 


5.7.2  Parameters 


Parameters  were  kept  at  their  default  values  with  the  following  exceptions : 


IV(MXFCAL)  -  min  {9999, 1000  *  n} 
IV(MXITER)  -  min  {9999, 1000  *  n} 
V(AFCTOL)  -  TOL2  (varied;  see  tables) 
V(RFCTOL)  -  TOL  (varied;  see  tables) 
V(  SCTOL)  -  eM 

V(  XCTOL)  -  TOL  (varied;  see  tables) 
V(  XFTOL)  -  eM 

V(  LMAXO)  -  usually  1.0  (default)  f 
V(  LMAXS)  -  1.0  (default) 

V(TUNERl)  -  0.1  (default) 


function  evaluation  limit 
iteration  limit 

absolute  function  convergence  tolerance 
relative  function  convergence  tolerance 
singular  convergence  tolerance 
x  convergence  tolerance 
false  convergence  tolerance 
initial  trust-region  diameter 
step  bound  for  singular  convergence  test 
reduction  test  coefficient 


f  In  some  cases  the  default  V(LMAXO)  =  1.0  for  the  initial  diameter  of  the  trust-region  was 
too  large  and  overflow  occurred  during  function  evaluation.  These  cases  are  indicated  in 
the  table  by  giving  the  lower  value  of  V(LMAXO)  that  was  subsequently  used  to  obtain  the 
results  in  the  column  labeled  “init.  diam.”. 

See  Dennis,  Gay,  and  Welsch  [1981a,  1981b],  Gay  [1983],  and  PORT  [1984]  for  details 
concerning  the  parameters. 


5.7.3  Convergence  Criteria 


The  following  quantities  will  be  used  in  describing  the  comvergence  criteria : 

objective  function  :  Tk  =  fjfk 
objective  gradient  :  gk  =  VTk  =  2 Jj  fk 

current  step  :  pk ,  the  minimizer  of  the  subproblem 

Newton  step  :  ^  «  P^tive  definite; 

l  undefined  otherwise. 

Newton  reduction  :  p.  =  (  if  ^  «  P°«ve  definite; 

I  0  otherwise, 

predicted  reduction  :  pF  =  —  Qk{Pk) 
actual  reduction  :  pA  =  Tk  —  T(xk  +  pk) 

..  -  V)Y\} 


Newton  reduction  :  pN 


_  |  -Qk{ps) 

1  0 


scaled  distance 


i/(x,y,D)  = 


maxi<t<„  {|(Dx)‘|  +  |(Dy)‘|}’ 


|  Here  v,  denotes  the  tth  component  of  the  vector  v.  There  is  a  provision  for  the 
user  to  replace  the  function  v\  we  used  the  default  in  all  of  the  tests. 

The  convergence  criteria  used  in  DN2G  are  as  follows : 

•  Absolute  function  convergence  occurs  at  x*  if 


ifki  <  v(afctol). 


(5.7.1) 


•  Relative  function  convergence  is  intended  to  approximate  the  condition 

Tk  ~  T(x,)  <  V(RFCTOL)  |^lfc| . 


The  test  actually  used  is 


pN  <  V(RFCTOL)  |jTfc| . 


(5.7.2) 


•  x  convergence  is  intended  to  approximate  the  condition 


v(xk,x\Dk )  <  V(XCTOL), 


The  test  actually  used  is 


Pk  =  Pn  and  v(xk,xk  +Pk,Dk)  <  V(XCTOL). 


(5.7.3) 


•  Singular  convergence  is  intended  to  approximate  the  condition 


Tk  -  min  { T(y )  |  \\Dk{y  -  a?*)||  <  v(uuxs)}  <  v(sctol)  \Tk\ , 
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where  Dk  is  the  diagonal  scaling  matrix  at  the  Jbth  iterate  —  when  none  of  the  con¬ 
vergence  criteria  listed  above  hold.  It  is  meant  to  indicate  relative  function  convegence 
when  the  Hessian  in  the  subproblem  is  singular. 

The  actual  test  is 

Tk  -  min  {Qfc(y)  |  ||Z)*(y  -  i*)||  <  v(LMAXS)}  <  v(sctol)  \Fk\ .  (5.7.4) 

Under  certain  conditions,  the  test  is  repeated  for  a  step  of  length  V(LMAXS). 

•  False  convergence  is  returned  if  none  of  the  other  convergence  criteria  are  satisfied 
and  a  trial  step  no  larger  than  V(XFCTOL)  is  rejected.  This  usually  indicates  either  an 
error  in  computing  the  objective  gradient,  or  a  discontinuity  (in  IF  or  g)  near  the  current 
iterate,  or  that  one  or  more  of  the  convergence  tolerances  (V(RFCTOL),  V(XCTOL),  and 
V(AFCTOL))  are  too  small  relative  the  accuracy  to  which  the  objective  is  computed. 
The  test  actually  used  is 

Fk  -  F(xk  +  Pk)  <  V(TUIEA1  )pP  and  xk  +  Pit, Dk)  <  v(XFTOL),  (5.7.5) 

where  the  parameter  V(TUNERl)  is  adjustable,  although  in  these  tests  the  default  value 
0.1  is  used  throughout. 

Except  for  test  (5.6.13),  tests  for  convergence  are  performed  only  when 

Pa  <  2 pP.  (5.7.6) 

See  Dennis,  Gay,  and  Welsch  [1981a,  1981b],  Gay  [1983],  and  PORT  [1984]  for  more 
discussion  of  the  convergence  criteria. 


The  following  abbreviations  are  used  in  the  tables  to  describe  the  conditions  under  which 
the  algorithm  terminates : 


ABS.  F  - 

(5.7.1) 

REL  F  - 

(5.7.2)  and  (5.7.6) 

X 

(5.7.3)  and  (5.7.6) 

X,  F  - 

(5.7.2)  and  (5.7.3)  and  (5.7.6) 

SING  - 

(5.7.4)  and  (5.7.6) 

(5.7.5)  and  (5.7.6) 

FALSE  - 

F  LIM  - 

function  evaluation  limit  reached 

TIME  - 

time  limit  exceeded 

LOOP  - 

subroutine  appears  to  loop 
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The  total  number  of  Jacobian  evaluations  is  either  equal  to  the  total  number  of  iterations 
of  the  method,  or  it  is  one  more  than  the  number  of  iterations.  The  number  in  the 
column  labeled  “iters.  /  J  evals.”  is  followed  by  a  “+"  if  an  extra  Jacobian  evaluation 
was  used  in  the  computation. 


90 


• 

♦ 

Numerical  Results  for  DN2G 

n  m  TOL 

init. 

/ 

iters./ 

ll*1l2 

ll/*ll2 

ll»*ll2 

est.  conv. 

diam. 

evals. 

J  evals. 

err. 

l.°  2  2  10"8 

14 

11+ 

1.41 

IQ-16 

IQ-15 

10-32  ABS.  r 

IO"12 

14 

11+ 

1.41 

IQ-16 

IQ-15 

IQ-32  ABS  p 

2.° 

2  2  10'8 

10 

8 

IO"12 

12 

1 

3.° 

2  2  IO-8 

IO'12 

4.° 

2  3  10“8 

40  I 

io-12 

53  4 

2  3 


10~8 

io-12 


6. 

2  10  IO"8 

14 

10+ 

io-12 

16 

12+ 

IO' 11 

o 

I 

o 

10-21 

ABS.  F 

0.00 

0.00 

0.00 

ABS.  F 

101 

10-2 

10"6 

REL.  ? 

101 

10-7 

IO-6 

REL.  P 

3  3 


IO-8 

IO"12 


9. 

3 

15 

10-8 

O 

1 

to 

10. 

CO 

16 

10-8  1 

10'12  1 

8 

IO"4 

io-12 

10-14 

8 

10"4 

Iq-16 

10“  14  X,  REL.  F 

3  10 


3  10 


4  4 


4  6 


4  11 


10“8 

10-12 


10-8 

io-12 


10“8 

IO'12 


10“8 

10-12 


10-8 

10-12 


32  120+ 

33  121 


(3000)  (2953) 

(3000)  (2953) 


X,  REL.  F 
X,  REL.  F 


16+ 

IO-4 

10'8 

o 

1 

»-* 

to 

io-17 

ABS.  P 

22+ 

10~6 

io-12 

1 

o 

10-24 

ABS.  P 

16. 

4 

20 

10"8 

10-12 

17. 

5 

33 

10-8 

IO"12 

18.° 

6 

13 

10'8 

10"12 

19. 

11 

65 

10'8 

10-12 

►— * 
O 

1 

<0 

i- 

1 

O 

10-17 

ABS.  F 

0.00 

0.00 

0.00 

ABS.  F 

10"2 

10-9 

10-9 

REL.  F 

IO'2 

10-9 

10-9 

REL.  P 

102 

IO"2 

10-8 

REL.  F 

102 

10-4 

IO-8 

REL.  F 

10-2 

IO"8 

10_n 

REL.  F 

IO"2 

10— !0 

10-11 

REL.  F 

10"9 

IO"8 

IO"17 

ABS.  F 

IO-16 

H- * 

0 

1 

<T> 

10-31 

ABS.  F 

IO"1 

IO-7 

10~8 

REL.  F 

IO"1 

IO'9 

10-8 

REL  P 

Numerical  Results  for  DN2G 


n 

m 

TOL 

init. 

/ 

iters./ 

ll**ll2 

II/II2 

n 

est. 

conv. 

diam. 

evals. 

J  evals. 

err. 

20a. 

6 

31 

10"8 

|| 

mm 

10"2 

10"9 

IO-10 

REL.  F 

lO"12 

K9 

BjW 

IO"2 

10'9 

IO-10 

REL.  F 

20b. 

9 

31 

10-8 

mm 

10-3 

10~u 

IO-13 

REL.  F 

io-12 

mm 

10-3 

10-14 

10-13 

X,  REL.  F 

20c. 

12 

31 

IO-8 

14 

n+ 

16.6 

IO"5 

10- 13 

10-ie 

REL.  F 

10-12 

14 

n+ 

16.6 

IO-5 

IO-13 

10"16 

REL.  F 

20d. 

20 

31 

10“8 

8 

6+ 

1.11 

10~8 

10-13 

10-16 

ABS.  F 

10-12 

(471) 

(137+) 

1.18 

10-8 

10-14 

10-16 

LOOP 

21a.° 

10 

10 

IO-8 

27 

17+ 

10-ie 

10-14 

10-si 

ABS.  F 

10-12 

27 

17+ 

10-ie 

IO-14 

10“31 

ABS-  F 

21b. 0 

20 

20 

10~8 

16 

12+ 

4.47 

10-16 

10-14 

10-31 

ABS.  F 

IO-12 

16 

12+ 

4.47 

10-16 

10-14 

10-si 

ABS.  F 

22a. 0 

12 

12 

10~8 

16+ 

IO"4 

10"8 

10-12 

IO-17 

ABS.  F 

10-12 

23+ 

IO-6 

10-12 

10-16 

10-24 

ABS.  F 

22b. 0 

20 

20 

IO-8 

20 

17+ 

IO"4 

10~8 

10-12 

IO-17 

ABS.  F 

IO-12 

26 

23+ 

10~6 

10~12 

10" 18 

10-24 

ABS.  F 

23a. 

4 

5 

10~8 

IO-3 

10-1° 

10-1° 

REL.  F 

10-12 

mSMl 

IO"3 

10-1° 

10-io 

REL.  P 

23b. 

10 

11 

IO"8 

61 

46+ 

10-2 

10-7 

IO-11 

REL.  F 

10-12 

68 

10-2 

10-10 

IO-11 

REL.  F 

24a. 

4 

8 

IO-8 

139 

IO-3 

IO"7 

IO-11 

REL.  P 

10-12 

142 

MEm 

10-3 

10"u 

10_n 

REL.  P 

24b. 

10 

20 

IO-8 

129 

mmm 

.598 

IO"2 

10-7 

10-9 

REL.  F 

IO-12 

138 

Kill 

.598 

IO-2 

10“8 

IO-9 

REL.  P 

25a. 0 

10 

12 

IO-8 

15 

10+ 

IO"12 

10"u 

10“24 

ABS.  F 

10-12 

16 

11+ 

10-15 

10-14 

10-31 

X 

25b. 0 

20 

22 

IO"8 

19 

12+ 

4.47 

IO"15 

IO"13 

IO"29 

X 

10-12 

19 

12+ 

4.47 

10-15 

10-13 

10“29 

ABS.  P 

26a. 0 

10 

10 

IO-8 

11 

7+ 

fennTiU 

IO-9 

IO-10 

IO-19 

ABS.  F 

10“12 

12 

8+ 

Wllil 

10-15 

10-i5 

10-30 

ABS.  F 

26b. 0 

20 

20 

IO"8 

10-3 

10-9 

10~6 

REL.  F 

IO-12 

IO-3 

10-1° 

IO-6 

REL.  F 

27a.° 

10 

10 

IO-8 

8 

6+ 

10-i0 

10-10 

IO-21 

ABS-  F 

IO-12 

9 

7+ 

10-15 

10“14 

IO-29 

ABS.  F 

27b. 0 

20 

20 

10“8 

mm 

m  | 

4.47 

IO-8 

10~8 

IO-17 

ABS.  F 

IO"12 

ma 

Mm 

4.47 

IO-14 

10-13 

IO-27 

ABS.  F 

28a.° 

10 

10 

10“8 

MM 

mm 

.412 

10-15 

10-16 

10-31 

ABS.  F 

10-12 

MM 

■■ 

.412 

10-15 

10-16 

10-31 

ABS.  P 

28b.° 

20 

20 

10~8 

mm 

IO-8 

10-9 

10“16 

ABS.  F 

IO-12 

10-16 

10-16 

IO-32 

ABS-  F 

92 


Numerical  Results  for  DN2G 


n 

m 

TOL 

init. 

diam. 

/ 

evals. 

iters./ 

J  ev&ls. 

nr  ii. 

IITII, 

eat. 

err. 

conv. 

29a. 0 

10 

10 

10'8 

mm 

mm 

.412 

10" 14 

IO-14 

10-29 

ABS.  P 

10-12 

mm 

KB 

.412 

10“14 

IO"14 

10- 29 

ABS.  P 

29b.° 

20 

20 

lO"8 

n 

mm 

10"14 

10~14 

10-28 

ABS.  P 

N 

w* 

1 

o 

•H 

H 

10- 14 

IO-14 

10-26 

ABS.  P 

30a.° 

10 

10 

IO-8 

8 

5+ 

IO-9 

IO"9 

10-16 

ABS.  P 

io-12 

9 

6+ 

2.05 

10- 16 

10-15 

IO"31 

ABS.  P 

30b. 0 

20 

20 

IO-8 

8 

5+ 

3.04 

10-8 

10-8 

►— * 
o 

1 

*4 

ABS.  P 

IO'12 

9 

6+ 

3.04 

10-15 

7 

o 

IQ-30 

ABS.  P 

31a. 0 

10 

10 

10"8 

10 

7+ 

1.80 

10-12 

IO"11 

10-23 

ABS.  P 

10-12 

11 

8+ 

1.80 

10-16 

10-15 

io-31 

X 

31b.° 

20 

20 

10"8 

10 

7+ 

2.66 

10-1° 

IO-10 

10-21 

ABS.  P 

10-12 

11 

8+ 

2.66 

10-15 

IO-14 

l0-3° 

ABS.  P 

32. L 

10 

20 

10~8 

5 

2 

3.16 

10° 

10-14 

IQ-16 

X,  REL.  P 

10-12 

5 

2 

3.16 

10° 

10~14 

IQ-16 

X,  REL.  P 

33. L 

10 

20 

IO"8 

18 

5 

20.2 

10° 

IO"11 

IO"6 

SINO. 

10-12 

18 

5 

20.2 

10° 

10_u 

10-6 

SINQ . 

34. L 

10 

20 

IO"8 

mm 

7.24 

10° 

10-1° 

10-6 

SING. 

10- 12 

mm 

7.24 

10° 

10-1° 

IO"6 

SING. 

35a. 

8 

8 

10-8 

mm 

1.65 

10"1 

10-7 

IO"9 

REL  P 

10-12 

■SB 

1.65 

10"1 

IO"8 

10~9 

REL.  P 

35b.° 

9 

9 

IO-8 

a 

7+ 

1.73 

10-12 

10-12 

10-24 

ABS.  P 

10- 12 

ii 

7+ 

1.73 

10-12 

IO-12 

10- 24 

ABS.  P 

35c. 

10 

10 

10-8 

17 

12+ 

1.81 

10" 1 

IO"7 

IO"3 

REL.  P 

IO-12 

19 

14+ 

1.81 

10"1 

IO"8 

IO-3 

REL.  P 

36a.° 

4 

4 

10“8 

(4000) 

(3988) 

27.7 

IO-7 

IO"7 

IO"13 

P  LIM. 

IO-12 

(4000) 

(3988) 

27.7 

10-7 

10-7 

10-13 

P  LIM. 

36b. 0 

9 

9 

IO"8 

(9000) 

(8977) 

35.5 

10-7 

10-8 

10- 14 

P  LIM. 

10-12 

(9000) 

(8977) 

35.5 

10-7 

IO"8 

10- 14 

P  LIM 

36c. 0 

9 

9 

IO"8 

16 

— 

1.73 

IO-9 

10-8 

10-17 

ABS.  P 

10-12 

22 

mSM 

1.73 

10-12 

10-12 

10-24 

ABS.  P 

36d.° 

9 

9 

IO-8 

(9000) 

(8966) 

38.4 

10-7 

10-7 

10" 14 

P  LIM. 

10-12 

(9000) 

(8966) 

38.4 

IO-7 

IO"7 

10-14 

P  LIM. 

37. 

2 

16 

IO-8 

10 

8+ 

8.85 

101 

IO"3 

IO-6 

REL.  P 

IO-12 

11 

9+ 

8.85 

101 

IO"5 

IO-6 

REL.  P 

38. 

3 

16 

IO"8 

10 

8+ 

26.1 

101 

10-2 

10-® 

REL.  P 

10-12 

12 

10+ 

26.1 

101 

10-7 

IO-6 

REL.  P 

93 


Numerical  Results  for  DH2G 


n 

m 

TOL 

init. 

diam. 

/ 

evals. 

iters./ 

J  evals. 

Ik*ll2 

lirila 

est. 

err. 

conv. 

39a. 

2 

3 

IO"8 

5 

4+ 

10"® 

IO"1 

10-1° 

IO-7 

REL-  P 

io-12 

5 

4+ 

IO"7 

10" 1 

IO-10 

IO"7 

REL.  P 

39b. 

2 

3 

IO-8 

kb 

IO-7 

IO-1 

10~7 

IO-7 

REL.  P 

io-12 

KB 

IO'7 

IO-1 

10~9 

IO"7 

REL.  F 

39c. 

2 

3 

IO"8 

7 

4+ 

10~7 

IO"1 

IO-7 

IO"7 

REL.  P 

IO"12 

8 

5+ 

IO-7 

IO-1 

10-1° 

10“7 

REL.  P 

39d. 

2 

3 

10~8 

7 

5+ 

IO"7 

IO-1 

10~® 

IO"7 

REL.  P 

IO-12 

8 

6+ 

10-7 

IO"1 

o 

1 

o 

IO"7 

REL.  P 

39e. 

2 

3 

IO-8 

9 

sa 

IO-8 

IO'1 

10-7 

IO"7 

REL.  P 

IO-12 

10 

mmi 

IO"8 

IO-1 

10-1° 

IO"7 

RBL.  P 

39f. 

2 

3 

IO-8 

14 

10+ 

IO-8 

IO'1 

10-6 

10“7 

REL.  P 

IO-12 

15 

11+ 

IO-9 

IO-1 

10~9 

10“7 

REL.  P 

39g. 

2 

3 

IO"8 

mmm 

IO"7 

IO-1 

10~4 

IO-7 

RBL.  P 

N 

1 

o 

rH 

Wm 

10~9 

IO"1 

10-® 

IO-7 

RBL.  P 

40a. 

3 

4 

10~8 

7 

6 

IO-6 

10° 

10-11 

IO"7 

RBL.  P 

10“12 

7 

6 

IO-6 

10° 

10- 11 

IO-7 

REL.  P 

40b. 

3 

4 

IO-8 

7 

5+ 

IO"6 

10° 

10-® 

10“7 

REL.  P 

10-12 

11 

9 

IO"6 

10° 

10~14 

IO-7 

REL.  P 

40c. 

3 

4 

IO"8 

9 

6+ 

10-7 

10° 

10-® 

IO-7 

REL.  P 

10-12 

10 

7+ 

IO-7 

10° 

10-9 

IO"7 

REL.  P 

40d. 

3 

4 

IO"8 

9 

7+ 

10-7 

10° 

10-8 

IO"7 

REL.  P 

10-12 

9 

7+ 

10-7 

10° 

10-8 

IO-7 

REL.  P 

40e. 

3 

4 

10“8 

10 

9+ 

10-7 

10° 

10-® 

IO-7 

REL  P 

10-12 

11 

10+ 

IO-7 

10° 

10“® 

IO"7 

RBL.  P 

40f. 

3 

4 

IO"8 

13 

10+ 

IO"8 

10° 

10-® 

IO"7 

REL.  P 

IO"12 

14 

11+ 

IO-8 

10° 

IO-7 

IO-7 

REL.  P 

40g. 

3 

4 

IO-8 

23 

16+ 

IO"8 

10° 

IO"4 

10-7 

REL.  P 

10-12 

25 

18+ 

IO"9 

10° 

IO"7 

IO"7 

REL.  P 

41a. 

5 

10 

IO-8 

4 

3+ 

IO"6 

10° 

10-10 

IO-7 

REL.  F 

10-12 

4 

3+ 

10-® 

10° 

10-1° 

IO-7 

REL.  P 

41b. 

5 

10 

10“8 

4 

3+ 

IO-6 

10° 

IO-7 

10-7 

REL.  P 

10-12 

5 

4+ 

IO-6 

10° 

10-11 

IO-7 

REL.  P 

41c. 

5 

10 

IO'8 

6 

5+ 

10“6 

10° 

10~8 

10~7 

REL.  P 

10-12 

6 

5+ 

10~® 

10° 

IO-8 

IO-7 

RBL.  P 

41d. 

5 

10 

IO-8 

9 

7+ 

10-® 

10° 

10-® 

IO-7 

RBL.  P 

IO-12 

11 

9+ 

10~® 

10° 

10~8 

IO-7 

REL.  P 

41e. 

5 

10 

IO-8 

13+ 

10“® 

10° 

10-® 

10“7 

REL.  P 

10- 12 

mm 

16+ 

IO"7 

10° 

IO-8 

IO-7 

REL.  P 

41f. 

5 

10 

10“8 

mm 

— 

10-® 

10° 

IO-4 

IO-7 

REL.  P 

IO-12 

mam 

IO"7 

10° 

IO-7 

IO-7 

REL.  P 

dig. 

5 

10 

IO-8 

IO-7 

10° 

10-5 

IO-7 

REL.  P 

10-12 

10~8 

10° 

10"® 

IO-7 

REL.  P 

94 


Numerical  Results  for  DH2G 


n 

m 

TOL 

init. 

di&m. 

/ 

e  vails. 

iters./ 

J  ev&ls. 

11**11, 

ll/*ll. 

Ha’ll, 

est. 

err. 

conv. 

42a. 0 

4 

24 

lO"8 

0.9 

29 

19+ 

60.8 

io-13 

10-1° 

10-25 

X 

io-12 

0.9 

29 

19+ 

60.8 

10-13 

10-1° 

IO-25 

ABS.  p 

42b. 0 

4 

24 

lO'8 

0.001 

74 

48+ 

61.3 

10-13 

10"a 

IO"25 

X 

10-12 

0.001 

74 

48+ 

61.3 

10-13 

10"n 

IO"25 

ABS.  P 

42c. 0 

4 

24 

IO-8 

0.01 

32 

19+ 

10-13 

10-1° 

IO"26 

X 

10"12 

0.01 

32 

19+ 

10-13 

10-1° 

10-28 

ABS.  P 

42d.° 

4 

24 

IO-8 

23 

18+ 

60.3 

10-10 

IO'7 

10“19 

ABS.  P 

10" 12 

24 

19+ 

60.3 

10- 13 

10-1° 

IO"26 

X 

43a.° 

5 

16 

IO-8 

10-12 

10-1° 

IO"24 

ABS.  P 

►— * 
o 

1 

10"14 

10_u 

IO"28 

X 

43b. 0 

5 

16 

10‘8 

20 

13+ 

54.0 

IO"12 

10-1° 

10"  24 

ABS.  P 

10-12 

20 

13+ 

54.0 

10-12 

10-1° 

10-24 

ABS.  P 

43c. 0 

5 

10-1 

10-7 

10-2 

REL.  P 

10"1 

10-1° 

10“2 

REL.  P 

43d. 0 

5 

16 

IO-8 

17 

n+ 

54.0 

10“14 

10-11 

10-27 

X 

10-12 

17 

n+ 

54.0 

10"14 

IO-11 

10-27 

ABS.  P 

43e.° 

5 

18+ 

mm i 

10"11 

10~9 

10-22 

ABS  P 

19+ 

10"14 

IO-12 

10-27 

X 

43f.° 

5 

16 

IO"8 

El 

15+ 

54.0 

10~12 

10- 10 

10-25 

ABS.  P 

10-12 

15+ 

54.0 

10-12 

10-1° 

10-25 

ABS.  P 

44a. 0 

6 

6 

10“8 

58 

41+ 

4.06 

10-1° 

TO 

1 

o 

IO-20 

ABS.  P 

w 

•H 

1 

o 

59 

42+ 

4.06 

10-15 

10-13 

10-3° 

ABS.  P 

44b. 0 

6 

6 

10'8 

■1 

SB 

10" 14 

10-12 

10- 28 

X 

10-12 

KB 

BB 

10"14 

10-12 

10-28 

ABS.  P 

44c. 0 

6 

6 

IO-8 

93 

84+ 

20.6 

10- 10 

10“6 

10-19 

ABS.  P 

IO"12 

94 

85+ 

20.6 

10-15 

IO"11 

10-29 

ABS.  P 

44d.° 

6 

6 

10~8 

97 

BE 

10"11 

10~8 

10-22 

ABS.  P 

10-12 

98 

IO-14 

IO"11 

10-28 

X 

44e.° 

6 

6 

10-8 

83 

72+ 

9.27 

10" 14 

10~12 

10-M 

X 

10-12 

83 

72+ 

9.27 

10"14 

10~12 

10-29 

ABS.  P 

45a. 0 

8 

8 

10~8 

65 

45+ 

4.06 

10-11 

10-9 

10-22 

ABS.  P 

10" 12 

66 

46+ 

4.06 

IO-17 

10-15 

10-33 

ABS.  P 

45b. 0 

8 

8 

10“8 

8 

7+ 

3.56 

IO*13 

IO'12 

10-26 

ABS.  P 

10-12 

8 

7+ 

3.56 

10-13 

10-12 

10-26 

ABS.  P 

45c. 0 

8 

8 

10~8 

129 

HE 

10-® 

10~4 

IO"16 

ABS.  P 

10-12 

130 

II 

10-15 

10-1° 

10-29 

ABS.  P 

45d.° 

8 

8 

IO-8 

168 

144+ 

15.3 

10- 14 

IO-11 

10-29 

X 

10-12 

168 

144+ 

15.3 

10-14 

IO-11 

10-29 

ABS.  P 

45e.° 

8 

8 

IO"8 

173 

165+ 

9.31 

10-15 

10-12 

10-29 

X 

10-12 

173 

165+ 

9.31 

10- 15 

10-12 

10-29 

ABS.  P 

95 


96 


5.8  Corrected  Gauss-Newton  Methods 

(NPL/NAG  LSQSDH  and  LSQFDQ) 


5.8.1  Software  and  Algorithms 

The  results  were  obtained  using  subroutines  LSQSDN  and  LSQFDQ  implementing  corrected 
Gauss- Newton  methods  from  the  National  Physical  Laboratory,  which  are  available  at  Stan¬ 
ford  Linear  Accelerator  Center.  A  subproblem  of  the  form 

min  gjp  +  +  Bk)p 

subject  toJtp  «  -  ft, 

is  solved  for  a  search  direction  pk,  where  ss  is  interpreted  in  a  least-squares  sense  using 
the  singular-value  decomposition  (see  Chapters  3  and  4).  Subroutine  LSQSDN  requires  exact 
second  derivatives  for  the  term  Bk  that  involves  the  second  derivatives  of  the  residuals, 
while  LSQFDQ  uses  a  quasi-Newton  approximation.  The  linesearch  algorithm  used  within  the 
subroutines  requires  both  function  and  gradient  information  (see  Gill  and  Murray  [1974], 
for  details).  These  subroutines  are  similar  to  those  available  in  the  NAG  Library  [1984]  for 
solving  nonlinear  least-squares  problems  :LSQSDN  corresponds  to  NAG  subroutine  E04HEF 
and  LSQFDQ  to  NAG  subroutine  E04GBF. 


5.8.2  Parmeters 


LSQSDN  and  LSQFDQ  have  the  same  set  of  input  parameters  as  the  corresponding  software 
from  the  NAG  Library  [1984],  The  values  chosen  are  listed  below. 


MAXCAL 

XTOL 

ETA 

STEPMX 


min  {9999, 1000  *  n} 
varied;  see  tables 
0.9 

usually  10s  (default)  f 


function  evaluation  limit 
accuracy  in  x 
linesearch  accuracy 
maximum  step  for  linesearch 


f  In  some  cases  the  default  STEPMX  =  10®  was  too  large  and  overflow  occurred  during  function 
evaluation  in  the  linesearch.  These  cases  are  indicated  in  the  tables  by  giving  the  lower  value  of 
STEPMX  that  was  subsequently  used  to  obtain  the  results. 

See  the  NAG  [1984]  manual  for  details  concerning  the  parameters. 
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5.8.3  Convergence  Criteria 


The  following  quantities  will  be  used  in  describing  the  convergence  criteria  : 


objective  function 
objective  gradient 
search  direction 
steplength 


n  =  flh 

9k  =  VTk  =  2  J?fk 

pk,  the  minimizer  of  the  subproblem 

ak,  determined  by  the  linesearch 


An  iterate  is  determined  to  be  optimal  by  LSQSDN  and  LSQFDQ  if 

<  cl 


(5.8.1) 


or 

llftlla  <  Cm||/*||2  (5.8.2) 

or  if  the  following  three  conditions  hold  : 

«*IIp*IIj  <  (MOL  +  c„)(l  +  ||*t||2)  (5.8.3) 

and 

F(*k-x)  -rk<  (XTOL  +  eM)2(l  +  |/*|)  (5.8.4) 

and 

M3<<l/3(1  +  I^l)-  (5-8.5) 

Conditions  (5.8.3)  and  (5.8.4)  are  meant  to  ensure  that  the  sequence  xk  has  converged, 
while  conditions  (5.8.2)  and  (5.8.5)  are  intended  to  test  whether  the  necessary  condition 
that  the  gradient  vanish  at  a  minimum  is  approximately  satisfied  at  xk.  Condition  (5.8.2) 
allows  the  algorithm  to  accept  a  point  as  a  local  mimimum  if  a  more  restrictive  test  on  the 
necessary  condition  than  (5.8.5)  is  met,  even  if  conditions  (5.8.3)  and  (5.8.4)  do  not  hold. 
For  the  zero-residual  case,  condition  (5.8.5)  specifies  that  the  method  may  also  terminate 
when  ||/t||2  is  no  larger  than  the  relative  machine  precision.  For  a  detailed  discussion  of 
convergence  criteria  similar  to  these,  see  Sections  8.2  and  8.5  of  Gill,  Murray,  and  Wright 
[1981].  In  particular,  Section  8.5. 1.3  treats  special  considerations  relevant  to  nonlinear  least 
squares. 
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The  following  abbreviations  are  used  in  the  tables  to  describe  the  conditions  under  which  the 
algorithm  terminates  : 


opt.  -  optimal  point  found 

*  -  current  point  cannot  be  improved  f 

F  lim.  -  function  evalutaion  limit  reached 
time  -  time  limit  exceeded 

t  A  corresponds  to  the  situation  in  which  the  algorithm  terminates  due  to  failure  in  the 
linesearch  to  find  an  acceptable  step  at  the  current  iteration. 
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Numerical  Results  for  LSQFDQ 
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10-9 

10"9 

IQ-16 

4c 

io-12 

11 

ma 

10"9 

10"9 

10"18 

4c 

31a.° 

10 

10 

10"8 

12 

8 

10"8 

10"7 

10-16 

4c 

10-12 

12 

8 

Bilil 

10"8 

10"7 

IQ-16 

* 

31b.° 

20 

20 

10~8 

12 

8 

2.66 

10"8 

10-7 

IO"16 

* 

IO"12 

12 

8 

2.66 

10"8 

10"7 

IO"16 

* 

32. 1 

10 

20 

10-8 

2 

1 

3.16 

10° 

10"14 

IO"16 

OPT. 

10-12 

2 

1 

3.16 

10° 

10~14 

IO"16 

OPT. 

33/ 

10 

20 

10~8 

2 

1 

■El 

10° 

10"9 

10-6 

OPT. 

10-12 

2 

1 

HB9 

10° 

10"9 

IO"6 

OPT. 

34/ 

10 

20 

10-8 

2 

1 

m a 

10° 

10"9 

IO"6 

OPT 

10-12 

2 

1 

10° 

10"9 

10"  6 

OPT. 

35a. 

8 

8 

10"8 

10.0 

mm 

mm 

10"1 

10"8 

IO"9 

4c 

10-12 

10.0 

KB 

Hi 

10"1 

10"8 

IO"9 

♦ 

35b.° 

9 

9 

10-8 

34 

14 

10-1° 

10-1° 

10-21 

* 

« 

1 

o 

•H 

34 

14 

IO"10 

10-1° 

IO"21 

♦ 

35c. 

10 

10 

10-8 

10.0 

mm 

10"1 

IO"9 

IO"3 

* 

10-12 

10.0 

an 

10"1 

10"9 

IO"3 

4c 

36a. 0 

4 

4 

10"8 

44 

25 

50.0 

10-15 

10-14 

10-51 

4c 

IO"12 

44 

25 

50.0 

10-15 

10"14 

10-31 

4c 

36b.° 

9 

9 

10"8 

44 

25 

10-15 

10"14 

IO"31 

4c 

10-12 

44 

25 

10-15 

10"14 

10-31 

* 

36c.° 

9 

9 

10-8 

28 

27 

10-15 

10" 16 

10-32 

OPT. 

IO"12 

28 

27 

10-16 

10-16 

IO"32 

OPT. 

36d.° 

9 

9 

10"8 

424 

100 

232. 

10-1° 

10"9 

10-20 

4c 

10-J2 

424 

100 

232. 

10-1° 

10"9 

10-2° 

* 

37. 

2 

16 

10-8 

25 

23 

8.85 

101 

10"4 

IO"6 

* 

. 

10-12 

25 

23 

8.85 

101 

10"4 

10'6 

4c 

38. 

3 

16 

10"8 

30 

26 

26.1 

101 

10"6 

IO"6 

+ 

o 

1 

>-* 

to 

30 

26 

26.1 

101 

10"6 

IO"6 

4c 

103 


Numerical  Results  for  LSQFDQ 


n  m  XTOL  max.  /,  J  iters.  ||x*||2  ll/'lla  ||j7*||2  est-  conv. 

step  evais.  err. 


39a. 

2 

3 

io-8 

io-12 

17 

17 

16 

16 

10“6 

IO-6 

*— *  M 

o  o 

1  1 

IO-9 

IO-9 

IO-7 

IO"7 

OPT. 

OPT. 

39b. 

2 

3 

IO-8 

24 

23 

10“7 

10" 1 

IO-9 

IO-7 

OPT. 

io-12 

24 

23 

IO"7 

10- 1 

10- 9 

IO-7 

OPT. 

39c. 

2 

3 

10~8 

22 

18 

IO-7 

IO-1 

IO"8 

IO-7 

OPT. 

IO'12 

23 

19 

10“7 

10- 1 

IO-9 

IO-7 

OPT. 

39d. 

2 

3 

10-8 

31 

26 

IO"7 

IO-1 

IO-8 

IO*7 

OPT. 

IO-12 

32 

27 

IO'7 

IO-1 

10~9 

IO"7 

OPT. 

39e. 

2 

3 

IO'8 

32 

26 

IO-8 

IO-1 

IO"8 

IO"7 

* 

10~12 

32 

26 

IO-8 

IO-1 

IO"8 

IO"7 

♦ 

39f. 

2 

3 

IO-8 

43 

27 

IO-9 

IO-1 

IO"8 

10-7 

* 

10-12 

43 

27 

IO-9 

10- 1 

IO-8 

IO-7 

* 

39g. 

2 

3 

10~8 

49 

28 

IO'9 

IO-1 

10~6 

IO-7 

* 

IO-12 

49 

28 

IO-9 

10- 1 

10~6 

IO'7 

* 

40a. 

3 

4 

IO"8 

18 

17 

10“6 

10° 

IO"8 

IO-7 

OPT. 

10-12 

18 

17 

IO-6 

10° 

10~8 

IO-7 

OPT. 

40b. 

3 

4 

IO"8 

19 

17 

IO-6 

10° 

IO-9 

IO-7 

OPT. 

10-12 

19 

17 

10~6 

10° 

IO"9 

IO"7 

OPT. 

40c. 

3 

4 

IO-8 

10~7 

10° 

IO-8 

IO-7 

OPT. 

IO-12 

10"  7 

10° 

nr8 

IO"7 

OPT. 

40d. 

3 

4 

10~8 

33 

26 

IO-7 

10° 

IO-8 

IO"7 

OPT. 

10- 12 

34 

27 

IO-7 

10° 

IO"9 

IO-7 

OPT. 

40e. 

3 

4 

IO"8 

— 

10-7 

10° 

IO-8 

10-7 

OPT. 

10-12 

mm 

IO"7 

10° 

IO"8 

IO-7 

♦ 

40f. 

3 

4 

IO-8 

92 

43 

IO-8 

10° 

IO-8 

IO"7 

* 

10-12 

92 

43 

IO-8 

10° 

IO"8 

IO"7 

* 

40g. 

3 

4 

IO"8 

IO3 

123 

53 

IO'9 

10° 

IO"6 

IO-7 

* 

10-12 

IO3 

123 

53 

IO"9 

10° 

IO-6 

IO'7 

* 

41a. 

5 

10 

10~8 

8 

7 

IO-6 

10° 

IO-9 

IO-7 

OPT. 

10-12 

8 

7 

IO-6 

10° 

IO-9 

IO-7 

OPT. 

41b. 

5 

10 

IO"8 

18 

13 

IO-6 

10° 

IO-9 

IO-7 

OPT. 

IO-12 

18 

13 

IO-6 

10° 

IO-9 

IO-7 

OPT. 

41c. 

5 

10 

IO"8 

21 

18 

IO-6 

10° 

IO-9 

IO"7 

OPT. 

10- 12 

21 

18 

10~6 

10° 

IO-9 

IO-7 

OPT. 

41d. 

5 

10 

IO-8 

38 

28 

10“6 

10° 

10“8 

IO-7 

* 

IO"12 

38 

28 

IO'6 

10° 

IO-8 

IO"7 

* 

41e. 

5 

10 

10~8 

47 

37 

IO-7 

10° 

IO-8 

IO"7 

* 

IO-12 

47 

37 

IO"7 

10° 

IO-8 

IO"7 

* 

41f. 

5 

10 

IO-8 

54 

45 

IO"7 

10° 

IO-8 

IO-7 

* 

10-12 

54 

45 

IO"7 

10° 

IO-8 

IO"7 

* 

41g. 

5 

10 

IO"8 

IO-8 

10° 

IO"7 

10~7 

★ 

10-12 

IO-8 

10° 

IO"7 

IO-7 

* 

104 


Numerical  Results  for  LSQFDQ 


n 

m 

XTOL 

max. 

step 

f,J 

evals. 

iters. 

m 

in 

est. 

err. 

conv. 

42a.° 

4 

24 

10"8 

10.0 

73 

49 

63.5 

IO"7 

10"“ 

IQ-13 

* 

10-12 

10.0 

73 

49 

63.5 

10~7 

IO"4 

10" 13 

* 

42b. 0 

4 

24 

10“8 

2.0 

94 

71 

61.9 

10-1° 

IO-7 

10-1® 

* 

IQ-12 

2.0 

94 

71 

61.9 

10-1° 

10“7 

10" 19 

* 

42c. 0 

4 

24 

10-8 

5.0 

46 

30 

60.3 

IO-5 

IO'3 

10“u 

* 

10"12 

5.0 

46 

30 

60.3 

IO'5 

IO-3 

IO'11 

♦ 

42d.° 

4 

24 

IO'8 

5.0 

60.3 

IO-5 

IO-3 

10-1° 

♦ 

io-12 

5.0 

60.3 

IO"5 

IO"3 

10-10 

♦ 

43a.° 

5 

16 

IO-8 

10.0 

33 

18 

54.0 

IO-9 

IO-6 

IO"17 

* 

io-12 

10.0 

33 

18 

54.0 

IO-9 

10“6 

10“17 

* 

43b. 0 

5 

16 

10'8 

10.0 

45 

21 

54.0 

IO-8 

IO"6 

IQ-16 

* 

10-12 

10.0 

45 

21 

54.0 

10~8 

10~6 

10- 16 

* 

43c.° 

5 

16 

IO”8 

10.0 

33 

18 

54.0 

IO-8 

10“6 

IO-16 

* 

IO-12 

10.0 

33 

18 

54.0 

IO-8 

IO"6 

10-18 

* 

43d.° 

5 

16 

IO-8 

10.0 

20 

54.0 

IO-9 

10“7 

IO*18 

4c 

10-12 

10.0 

20 

54.0 

IO"9 

IO"7 

10-18 

* 

43e.° 

5 

16 

IO-8 

10.0 

27 

14 

54.0 

l0-i° 

IO"8 

10-2° 

* 

IO-12 

10.0 

27 

14 

54.0 

10-ic 

IO-8 

10- 20* 

* 

43f.° 

5 

16 

IO"8 

mm 

i  a  i 

54.0 

10~8 

IO"6 

10-1« 

♦ 

10-12 

■  ■  ■ 

K9 

54.0 

10~8 

IO"6 

10- 1* 

♦ 

44a.° 

6 

6 

IO-8 

97 

28 

4.06 

io-12 

IO-10 

IO"23 

4C 

IO-12 

97 

28 

4.06 

10-12 

10- 10 

IO-23 

4c 

44b. 0 

6 

6 

10~8 

10 

6 

3.52 

IO"8 

10-6 

10-15 

4c 

10~12 

10 

6 

3.52 

IO"8 

10“6 

10-15 

* 

44c. 0 

6 

6 

IO'8 

47 

19 

20.6 

IO"7 

10-3 

10-14 

* 

10-12 

47 

19 

20.6 

10-7 

10"3 

io-14 

♦ 

44d.° 

6 

6 

IO"8 

40 

17 

15.3 

10-9 

IO"5 

IO-17 

4c 

IO"12 

40 

17 

15.3 

10"9 

IO-5 

IO-17 

* 

44e.° 

6 

6 

IO*8 

47 

20 

9.27 

10- 13 

10-10 

IO-25 

4c 

10-17 

47 

20 

9.27 

10“ 13 

10-1° 

10-25 

4c 

45a.° 

8 

8 

IO"8 

97 

28 

4.06 

IO"12 

IO-10 

IO-23 

* 

10-12 

97 

28 

10-12 

10-1° 

10~23 

4c 

45b. 0 

8 

8 

IO'8 

12 

mm 

IO-8 

IO-6 

10-15 

4c 

10-12 

12 

mm 

10“8 

IO-6 

10-15 

4c 

45c.° 

8 

8 

IO-8 

47 

19 

20.6 

IO"7 

IO-3 

10~14 

* 

10-12 

47 

19 

20.6 

10-7 

IO-3 

IO-14 

4c 

45d.° 

8 

8 

IO-8 

42 

18 

15.3 

IO-9 

IO-5 

IO-17 

4c 

10-12 

42 

18 

15.3 

10“9 

IO"5 

IO-17 

* 

45e.° 

8 

8 

IO-8 

49 

21 

9.31 

10-13 

10-1° 

10-25 

4c 

IO-12 

49 

21 

9.31 

IO-13 

10-1° 

10-25 

* 

105 


Numerical  Results  for  LSQSDV 


n 

m 

XTOL 

max. 

step 

f,J 

e  veils. 

iters. 

11**11, 

ll/*ll. 

lls*ll. 

est. 

err. 

conv. 

l.° 

2 

2 

10-8 

31 

12 

1.41 

0.00 

0.00 

0.00 

OPT. 

io-1J 

31 

12 

1.41 

0.00 

0.00 

0.00 

OPT. 

2.° 

2 

2 

10"8 

16 

6 

11.4 

101 

10-9 

101 

OPT. 

10"12 

18 

8 

11.4 

101 

10-® 

101 

OPT. 

3.° 

2 

2 

10"8 

47 

25 

9.11 

10-10 

10- 5 

10-20 

* 

10-12 

47 

25 

9.11 

10“ 10 

10"® 

10"  20 

* 

4.° 

2 

3 

10-8 

53 

21 

10* 

10-9 

10-3 

10-ia 

* 

10-12 

53 

21 

10* 

10"9 

10-3 

10-18 

* 

5.° 

2 

3 

10"8 

mm 

mm 

10" 14 

10-13 

10-28 

OPT. 

10-12 

■a 

SB 

10" 14 

10“ 13 

10-2S 

OPT. 

6. 

2 

10 

10"8 

5.0 

36 

10 

.365 

101 

10“5 

10-6 

♦ 

10-12 

5.0 

36 

10 

.365 

101 

10“5 

10-® 

* 

7.° 

3 

3 

10“8 

14 

10 

1.00 

10-11 

10-10 

10-28 

OPT. 

10-12 

14 

10 

1.00 

10"11 

10- 10 

1Q-2S 

OPT. 

8. 

3 

15 

10-8 

6 

5 

2.60 

10"1 

10~9 

10-8 

OPT. 

10-12 

6 

5 

2.60 

10-1 

10-9 

10-8 

OPT. 

9. 

3 

15 

10"8 

3 

2 

1.08 

10"4 

10-12 

10"M 

OPT. 

>-* 

o 

1 

to 

3 

2 

1.08 

10-4 

10-12 

10~14 

OPT. 

10. 

3 

16 

10~8 

17 

10 

104 

101 

102 

10"® 

* 

10- 12 

17 

10 

104 

10* 

102 

10-* 

* 

11.° 

3 

10 

10“8 

10.0 

30 

15 

10-2 

10-9 

10-8 

* 

10-12 

10.0 

30 

15 

KWH 

10"2 

10“9 

10-8 

* 

12.° 

3 

10 

10"8 

8 

6 

10.1 

10-1° 

10- 10 

10-19 

OPT. 

10-12 

12 

8 

10.1 

10-1° 

10-10 

10-19 

♦ 

13.° 

4 

4 

10-8 

18 

17 

10-5 

10-9 

10-is 

10-is 

OPT. 

10-12 

18 

17 

10"5 

10"9 

10-13 

10-18 

OPT. 

14.° 

4 

6 

10-8 

87 

42 

mm 

10-8 

10-7 

10-16 

OPT 

10-12 

93 

45 

1 

10'8 

10-7 

10-16 

* 

15. 

4 

11 

10-8 

16 

8 

.328 

10-2 

10- 14 

10~9 

OPT. 

K- » 

o 

1 

to 

16 

8 

.328 

10-2 

10-14 

10"9 

OPT. 

16. 

4 

20 

10"8 

45 

22 

17.6 

102 

10-8 

10-8 

* 

10- 12 

45 

22 

17.6 

102 

10"® 

10"8 

* 

17. 

5 

33 

10-8 

■9 

9 

mm 

10"2 

10-9 

10- 11 

OPT. 

10-12 

mm 

11 

EH 

10*2 

10-9 

10-u 

* 

O 

00 

r-4 

6 

13 

10-8 

10.0 

243 

69 

12.3 

10-12 

10“u 

10-23 

OPT. 

10-12 

10.0 

247 

71 

12.3 

10-12 

10-" 

10-23 

♦ 

19. 

11 

65 

10-8 

19 

10 

9.38 

10-1 

10“9 

10-8 

OPT. 

10"12 

19 

10 

9.38 

10-* 

10"9 

10-8 

OPT. 
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Numerical  Results  for  LSQSDH 


n 

m 

XTOL 

max. 

step 

f,J 

evals. 

iters. 

11**11* 

urn. 

n 

est. 

err. 

conv. 

20a. 

6 

31 

ur8 

9 

7 

2.44 

io-2 

IO"8 

10-10 

* 

10" 12 

9 

7 

2.44 

io-2 

10~8 

IO-10 

* 

20b. 

9 

31 

10-8 

6 

5 

6.06 

10"3 

IO-11 

IQ-13 

OPT. 

lO"12 

6 

5 

6.06 

10-3 

10“n 

10-18 

OPT. 

20c. 

12 

31 

IO"8 

6 

5 

16.6 

IO"5 

10~u 

10-16 

OPT. 

o 

1 

►— 

to 

6 

5 

16.6 

IO"5 

10“ 14 

10-16 

OPT. 

20d. 

20 

31 

KT8 

10 

8 

247. 

10-1° 

10"n 

10-24 

OPT. 

lO"12 

12 

9 

247. 

IO"10 

IO"11 

10-24 

OPT. 

21a. 0 

10 

10 

lO-8 

31 

12 

3.16 

0.00 

0.00 

0.00 

OPT. 

10"12 

31 

12 

3.16 

0.00 

0.00 

0.00 

OPT. 

21b.° 

20 

20 

10“8 

31 

12 

4.47 

0.00 

0.00 

0.00 

OPT. 

io-12 

31 

12 

4.47 

0.00 

0.00 

0.00 

OPT. 

22a. 0 

12 

12 

10~8 

18 

17 

IO-5 

10~9 

IO"13 

10-is 

OPT. 

10-12 

18 

17 

IO"5 

IO-9 

IO"13 

10-18 

OPT. 

22b.° 

20 

20 

10"8 

18 

17 

IO-5 

IO"9 

IO"13 

IO-18 

OPT. 

10-12 

18 

17 

10"  5 

IO"9 

10-13 

10-18 

OPT. 

23a. 

4 

5 

10~8 

58 

32 

.500 

IO-3 

10-“ 

10-1° 

OPT. 

>— » 
o 

1 

58 

32 

.500 

10~3 

IO-11 

IO-10 

OPT. 

23b. 

10 

11 

10"8 

124 

64 

.500 

IO"2 

10- 13 

10" 11 

OPT. 

10-12 

124 

64 

.500 

10-2 

10-13 

10- 11 

OPT. 

24a. 

4 

8 

10"8 

228 

136 

.759 

10-3 

10-11 

10-" 

OPT. 

k— * 

o 

1 

►-* 

10 

228 

136 

.759 

IO”3 

10-11 

10-11 

OPT. 

24b. 

10 

20 

10'8 

88 

.598 

10-2 

IO"11 

IO"9 

OPT. 

k—» 

o 

1 

88 

.598 

IO"2 

10-11 

10-9 

OPT. 

25a.° 

10 

12 

1(T8 

12 

10 

3.16 

IO"8 

10“® 

10-15 

OPT. 

10"12 

17 

12 

3.16 

IO"8 

10~6 

10-15 

♦ 

25b.° 

20 

22 

10"8 

14 

12 

4.47 

IO"8 

IO"7 

10“17 

OPT. 

10-12 

19 

14 

4.47 

IO"8 

IO"7 

IO-17 

* 

26a.° 

10 

10 

1 

o 

•-4 

18 

9 

.306 

10_u 

IO-11 

10-22 

OPT. 

IO"12 

22 

11 

.306 

10"n 

10"u 

10-22 

* 

26b.° 

20 

20 

10“8 

18 

9 

.208 

10-10 

10-1° 

IO-19 

OPT. 

10-12 

26 

12 

.208 

10-1° 

10-1° 

10-1* 

* 

27a.° 

10 

10 

IO-8 

102 

22 

7 

3.18 

10-1° 

IO"9 

10-2° 

OPT. 

10-12 

102 

28 

10 

3.18 

10-1° 

10~9 

10-20 

* 

27b.° 

20 

20 

10"8 

102 

21 

7 

4.47 

IO-9 

IO"9 

IO"19 

OPT. 

10-12 

102 

27 

10 

4.47 

IO-9 

IO"9 

IO-19 

* 

28a. 0 

10 

10 

IO-8 

mm 

.412 

10-15 

10-16 

10-81 

OPT. 

fl 

1 

© 

.412 

10-15 

10-16 

10-31 

OPT. 

28b.° 

20 

20 

10"8 

n 

.571 

10-16 

10-16 

10-32 

OPT. 

10-12 

WM 

.571 

10-16 

10-16 

IO- 32 

OPT. 
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Numerical  Results  for  LSQSDtf 


n 

m 

XTOL 

m&x. 

step 

f,J 

evals. 

iters. 

11**11, 

ll/*ll, 

11*11, 

est. 

err. 

conv. 

29a.° 

10 

10 

IO-8 

6 

4 

.412 

io-14 

IO"14 

10-2® 

OPT. 

ur12 

6 

4 

.412 

10-14 

10~14 

IO"29 

OPT. 

29b.° 

20 

20 

IO"8 

6 

4 

.571 

io-14 

10“14 

10-2S 

OPT. 

IO"12 

6 

4 

.571 

10"14 

10-14 

10-26 

OPT. 

30a. 0 

10 

10 

IO-8 

7 

IO-9 

IO-9 

10-16 

OPT. 

1 

o 

rH 

11 

10-9 

IO-9 

IO-!8 

♦ 

30b.° 

20 

20 

10"8 

7 

5 

3.04 

IO"9 

IO"9 

IQ-16 

OPT. 

IO"12 

11 

7 

3.04 

IO-9 

10‘9 

IO'18 

* 

31a.° 

10 

10 

IO-8 

8 

6 

1.80 

IO-8 

10“7 

IQ-16 

OPT. 

10-12 

12 

8 

1.80 

10~8 

10-7 

10-16 

♦ 

31b.° 

20 

20 

10"8 

8 

6 

2.66 

IO"8 

IO-7 

IQ-16 

OPT. 

IO"12 

12 

8 

2.66 

IO-8 

IO"7 

IQ-16 

* 

32.L 

10 

20 

10-8 

2 

1 

3.16 

10° 

10-14 

IQ-16 

OPT. 

10-12 

2 

1 

3.16 

10° 

IO-14 

IQ-16 

OPT. 

33. 1 

10 

20 

10-8 

2 

1 

BE9 

10° 

10-9 

10~6 

OPT. 

IO-'2 

2 

1 

m 

10° 

IO"9 

10~6 

OPT. 

34. 1 

10 

20 

IO-8 

2 

1 

Bi 

10° 

10-9 

10~6 

OPT 

10-12 

2 

1 

mSm 

10° 

10"® 

10~6 

OPT. 

35r. 

8 

8 

no 

1 

o 

74 

19 

1.65 

10-1 

IO"7 

10~9 

* 

o 

1 

to 

74 

19 

1.65 

10- 1 

IO"7 

10~9 

♦ 

35b.° 

9 

9 

10"8 

30 

12 

10-“ 

IO-10 

IO"22 

OPT. 

10- 12 

34 

14 

IO-11 

IO-10 

10-22 

* 

35c. 

10 

10 

10-8 

43 

10 

1.81 

10-1 

10-10 

10~3 

OPT. 

10-12 

43 

10 

1.81 

10-1 

IO"10 

10~3 

OPT. 

36a.° 

4 

4 

10"8 

38 

22 

50.0 

10-15 

IO-14 

IO"31 

OPT. 

10-12 

38 

22 

50.0 

10-15 

10-14 

IO-31 

OPT. 

36b.° 

9 

9 

IO-8 

38 

22 

50.0 

10-15 

IO-14 

IO"31 

OPT. 

10-12 

38 

22 

50.0 

10-15 

IO'14 

IO"31 

OPT. 

27 

10-16 

10-16 

IO-32 

OPT. 

27 

10-16 

10-16 

IO-32 

OPT. 

36d.° 

9 

9 

IO-8 

ES9 

84 

233. 

10-10 

IO"9 

IO-20 

OPT. 

10-12 

EH 

84 

233. 

10-10 

IO"9 

IO-20 

OPT. 

37. 

2 

16 

10~8 

9 

8 

8.85 

101 

10-1° 

10~6 

OPT. 

10-12 

9 

8 

8.85 

101 

o 

1 

o 

»*H 

10-6 

OPT. 

38. 

3 

16 

IO-8 

10 

8 

26.1 

101 

IO-9 

10~6 

OPT. 

10-12 

10 

8 

26.1 

101 

10-9 

10~6 

OPT. 
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Numerical  Results  for  LSQSDN 


n 

m 

XTOL 

max. 

step 

f,J 

evals. 

iters. 

n 

in 

ll*‘ll2 

est. 

err. 

conv. 

39a. 

2 

3 

lO"8 

4 

3 

IO*6 

IO"1 

IO'9 

10~7 

OPT. 

i<r12 

4 

3 

1(T® 

IO-1 

IO-9 

IO"7 

OPT. 

39b. 

2 

3 

IO-8 

6 

5 

IO"7 

IO"1 

IO-15 

IO"7 

OPT. 

10-12 

6 

5 

IO"7 

IO-1 

10-!5 

IO-7 

OPT. 

39c. 

2 

3 

IO"8 

9 

5 

10-7 

IO'1 

10-1° 

KT7 

OPT. 

IO-12 

9 

5 

10-7 

IO-1 

10-10 

IO-7 

OPT. 

39d. 

2 

3 

IO"8 

12 

7 

10-7 

IO-1 

10-10 

IO-7 

OPT. 

IO-12 

12 

7 

10-7 

IO-1 

10-10 

IO-7 

OPT. 

39e. 

2 

3 

10“8 

12 

7 

IO-8 

IO-1 

IO-10 

IO"7 

OPT. 

IO-12 

12 

7 

IO"8 

IO-1 

IO"10 

10“7 

OPT. 

39f. 

2 

3 

IO-8 

25 

10 

IO"9 

IO-1 

IO-13 

10“7 

OPT. 

IO-12 

25 

10 

IO"9 

IO-1 

10-13 

IO-7 

OPT. 

39g. 

2 

3 

IO-8 

mm 

10-1° 

IO-1 

IO"7 

IO-7 

* 

IO-12 

89 

10-10 

IO-1 

IO-7 

10- 7 

* 

40a. 

3 

4 

1(T8 

5 

4 

10-® 

10° 

IO-15 

IO-7 

OPT. 

IO-12 

5 

4 

IO"6 

10° 

IO"15 

IO"7 

OPT. 

40b. 

3 

4 

6 

4 

IO-6 

10° 

10-10 

IO-7 

OPT. 

IO-12 

6 

4 

IO"6 

10° 

10-!O 

IO"7 

OPT. 

40c. 

3 

4 

11 

6 

IO"7 

10° 

IO"14 

IO-7 

OPT. 

IO"12 

11 

6 

IO"7 

10° 

IO"14 

IO"7 

OPT. 

40d. 

3 

4 

13 

7 

IO"7 

10° 

10- 16 

IO"7 

OPT. 

10- 12 

13 

7 

IO-7 

10° 

10-16 

IO"7 

OPT. 

40e. 

3 

4 

45 

16 

10-7 

10° 

IO-13 

IO-7 

OPT. 

IO-12 

45 

16 

IO"7 

10° 

10~13 

IO-7 

OPT. 

40f. 

3 

4 

49 

18 

10~8 

10° 

10-15 

IO-7 

OPT. 

IO-12 

49 

18 

IO-8 

10° 

10-15 

IO-7 

OPT. 

40g. 

3 

D 

a 

IO"9 

10° 

t— * 

o 

1 

o 

IO-7 

OPT. 

■ 

19 

IO-9 

10° 

10-1° 

IO-7 

OPT. 

41a. 

5 

4 

3 

IO-6 

10° 

10-12 

IO-7 

OPT. 

IO"12 

4 

3 

10“8 

10° 

10-12 

IO"7 

OPT. 

41b. 

5 

4 

3 

IO-6 

10° 

IO"9 

IO-7 

OPT. 

10-n 

4 

3 

IO-6 

10° 

IO-9 

IO"7 

OPT. 

41c. 

5 

10 

IO-8 

5 

4 

10"6 

10° 

10-9 

IO"7 

OPT. 

10- 12 

5 

4 

IO"6 

10° 

IO-9 

IO-7 

OPT. 

41d. 

5 

9 

7 

10-® 

10° 

IO"12 

IO-7 

OPT. 

IO-12 

9 

7 

10”® 

10° 

10-12 

IO-7 

OPT. 

41e. 

5 

10 

IO-8 

14 

10 

IO-7 

10° 

IO-11 

IO'7 

OPT. 

10-12 

14 

10 

IO-7 

10° 

10- 11 

IO-7 

OPT. 

41f. 

5 

10 

IO"8 

16 

12 

10-7 

10° 

10-1° 

IO-7 

OPT. 

IO"12 

16 

12 

IO"7 

10° 

10-10 

IO-7 

OPT. 

41g. 

5 

10 

IO-8 

mm 

mm 

IO-8 

10° 

IO-12 

IO-7 

OPT- 

IO"12 

89 

89 

IO-8 

10° 

IO"12 

IO"7 

OPT. 
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Numerical  Results  for  LSQSDV 


n 

m 

XTOL 

max. 

step 

f.J 

evals. 

iters. 

urn. 

est. 

err. 

conv. 

42a. 0 

4 

24 

10'8 

mm 

mm 

63.5 

lO-io 

10-7 

10-i® 

OPT. 

10"12 

■9 

mm 

63.5 

lO-io 

10-7 

10-19 

♦ 

42b. 0 

4 

24 

10“8 

2.0 

75 

67 

61.9 

10~6 

10-3 

10-11 

* 

10-12 

2.0 

75 

67 

61.9 

10~6 

10“3 

10-n 

* 

42c. 0 

4 

24 

10"8 

mm 

48 

28 

61.8 

10-12 

10-® 

10-23 

OPT. 

10-12 

HI 

48 

28 

61.8 

10-12 

10- 9 

10-23 

OPT. 

42d.° 

4 

24 

10-8 

5.0 

27 

19 

mmm 

10~5 

10"3 

10"10 

* 

10"12 

5.0 

27 

19 

10~5 

10"3 

10-1° 

* 

43a. 0 

5 

16 

10~8 

10.0 

25 

14 

54.0 

10~® 

10-® 

10-17 

OPT 

C4 

H 

1 

o 

10.0 

33 

18 

54.0 

10~9 

10"6 

10~17 

* 

43b. 0 

5 

16 

lO-8 

10.0 

37 

17 

54.0 

10~8 

10-® 

10-1® 

OPT. 

o 

1 

h* 

NJ 

10.0 

45 

21 

54.0 

10~8 

10-® 

10-1® 

* 

43c. 0 

5 

16 

10~8 

10.0 

25 

14 

54.0 

10~8 

10-® 

10-16 

OPT. 

10'12 

10.0 

33 

18 

54.0 

10~8 

10"® 

10“16 

* 

43d. 0 

5 

16 

10"8 

10.0 

16 

54.0 

10~9 

10~7 

10" 18 

OPT. 

10-12 

10.0 

38 

54.0 

10~® 

10-7 

10- 18 

* 

43e.° 

5 

16 

10"8 

10.0 

19 

10 

54.0 

10-1° 

10"8 

10-2° 

OPT. 

10-12 

10.0 

27 

14 

10-1° 

10-8 

10-20 

* 

43f.° 

5 

16 

10-8 

cm 

10~8 

10-® 

10-i® 

OPT. 

10-12 

■9 

EM 

10~8 

10"® 

10" 16 

* 

44a.° 

6 

6 

10“8 

93 

26 

10- 12 

10-10 

10-23 

OPT. 

10-12 

95 

27 

4.06 

10-12 

10-10 

10-23 

OPT. 

44b. 0 

6 

6 

10“8 

6 

4 

3.52 

10'8 

10"® 

10-15 

OPT. 

10-12 

10 

6 

3.52 

10~8 

10"® 

10-15 

* 

44c. 0 

6 

6 

10-8 

mm 

mm 

20.6 

10-7 

10-3 

10“14 

* 

10-12 

H 

KB 

20.6 

10-7 

10-3 

10-14 

* 

44d.° 

6 

6 

10~8 

40 

17 

15.3 

10~® 

10"5 

10-17 

* 

10-12 

40 

17 

15.3 

10~® 

10-® 

10"17 

* 

44e.° 

6 

6 

10~8 

BE 

10-13 

10-1° 

10-25 

OPT. 

10-12 

SI 

10"13 

10-10 

10"25 

OPT. 

45a. 0 

8 

8 

10"8 

93 

26 

Bia 

10-12 

10-10 

10-23 

OPT. 

10-12 

95 

27 

RH 

10-12 

10-10 

10-23 

OPT. 

45b. 0 

8 

8 

10-8 

6 

4 

3.56 

10~8 

10-® 

10-15 

OPT. 

10"12 

12 

7 

3.56 

10~8 

10-® 

10-15 

* 

45c. 0 

8 

8 

10"8 

47 

10-7 

10"3 

10"14 

♦ 

10-12 

47 

10~7 

10-3 

10~14 

* 

45d.° 

8 

8 

10-8 

42 

18 

15.3 

10~® 

10-® 

10-17 

* 

<N 

H 

1 

o 

42 

18 

15.3 

10~9 

10-® 

10-17 

* 

45e.° 

8 

8 

10-8 

mm 

9.31 

10- 13 

10"10 

10-2S 

OPT 

10-12 

mlM 

9.31 

10- 13 

10-1° 

10-25 

OPT. 
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5.9  Test  Problems 


♦ 


Superscripts  on  problem  numbers  have  the  following  interpretation  : 

0  :  zero-residual  problem 
1  :  linear  least-squares  problem 


Problems  from  More,  Garbow,  and  Hillstrom  [1981] 


l.° 

n 

2 

m 

2 

Rosenbrock 

2.° 

2 

2 

Freudenstein  and  Roth 

3.° 

2 

2 

Powell  Badly  Scaled 

4.° 

2 

3 

Brown  Badly  Scaled 

5.° 

2 

3 

Beale 

6. 

2 

10 

Jennrich  and  Sampson 

7.° 

3 

3 

Helical  Valley 

8. 

3 

15 

Bard 

9. 

3 

15 

Gaussian 

10. 

3 

16 

Meyer 

11.° 

3 

10 

Gulf  Research  and  Development] 

12.° 

3 

10 

Box  3-Dimensional 

13.° 

4 

4 

Powell  Singular 

14.° 

4 

6 

Wood 

15. 

4 

11 

Kowalik  and  Osborne 

16. 

4 

20 

Brown  and  Dennis 

17. 

5 

33 

Osborne  1 

18.° 

6 

13 

Biggs  EXP6J 

f  For  the  Gulf  Research  and  Development  Function  (#  11),  the  formula 

f  \y*mix2\x 3‘ 


<t>i(x)  =  exp 


x i 


-U 


given  in  More,  Garbow,  and  Hillstrom  [1981]  for  the  residual  functions  is  in  error.  The 
correct  formula  is 

<t>i{x)  =  exp 

(see  More,  Garbow,  and  Hillstrom  [1978]). 


lift  ~  »2l 


*3 


l! 


J  For  the  Biggs  EXP6  Function  (#  18),  the  minmum  value  for  the  sum  of  squares 
is  given  in  Mori,  Garbow,  and  Hillstrom  [1981]  as  5.65565.. .  X  10-3.  It  can  be  easily 
verified  that  the  residuals  vanish  at  several  points  (for  example  (1,10,1,5,4,3)). 
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Problems  from  More,  Garbow,  and  Hillstrom  [1981]  (continued) 


19. 

n 

11 

m 

65 

Osborne  2f 

20a. 

6 

31 

Watson 

20b. 

9 

31 

Watson 

20c. 

12 

31 

Watson 

20d. 

20 

31 

Watson 

21a.° 

10 

10 

Extended  Rosenbrock 

21b.° 

20 

20 

Extended  Rosenbrock 

22a.° 

12 

12 

Extended  Powell  Singular 

22b.° 

20 

20 

Extended  Powell  Singular 

23a. 

4 

5 

Penalty  I 

23b. 

10 

11 

Penalty  I 

24a. 

4 

8 

Penalty  II 

24b. 

10 

20 

Penalty  II 

25a.° 

10 

12 

Variably  Dimensioned 

25b.° 

20 

22 

Variably  Dimensioned 

26a.° 

10 

10 

Trigonometric 

26b.° 

20 

20 

Trigonometric 

27a.° 

10 

10 

Brown  Almost  Linear 

27b.° 

20 

20 

Brown  Almost  Linear 

28a.° 

10 

10 

Discrete  Boundary  Value 

28b.° 

20 

20 

Discrete  Boundary  Value 

29a.° 

10 

10 

Discrete  Integral 

29b.° 

20 

20 

Discrete  Integral 

30a.° 

10 

10 

Broyden  Tridiagonal 

30b.° 

20 

20 

Broyden  Tridiagonal 

31a.° 

10 

10 

Broyden  Banded 

31b.° 

20 

20 

Broyden  Banded 

32. 1 

10 

20 

Linear — Full  Rank 

33. 1 

10 

20 

Linear — Rank  1 

34/ 

10 

20 

Linear — Rank  1  with  Zero  Rows  and  Cols, 

35a. 

8 

8 

Chebyquad 

35b.° 

9 

9 

Chebyquad 

35c. 

10 

10 

Chebyquad 

f  For  Osborne’s  Second  Function  (#  19),  the  value  of  /(**)  is  given  (to  six  figures) 
in  Mor£,  Garbow,  and  Hillstrom  [1981]  as  4.01377  X  10-2.  The  smallest  value  we  were 
able  to  obtain  was  4.01683  X  10-2. 
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Matrix  Square  Root  Problemsft 


n  m 
36a.°  4  4 

36b.°  9  9 

36c.°  9  9 

36d.°  9  9 


Matrix  Square  Root  1 
Matrix  Square  Root  2 
Matrix  Square  Root  3 
Matrix  Square  Root  4 


f  These  test  problems  come  from  a  private  communication  of  S.  Hammarling  to  P.  E. 
Gill  in  1983. 


MATRIX 


SQUARE  ROOT 


36a. 


36b. 


36c. 


36d. 


fl 0~4  1  \  /10"2  50  ^ 

V  o  io~4)  V  0  10-2  / 


/HT4 

1 

0  \ 

/10"2 

50 

0  \ 

(  0 

10"4 

0  I 

I  0 

io-2 

0  ) 

V  o 

0 

10 “V 

\  0 

0 

10-2/ 

}  The  identity  matrix  was  used  as  the  starting  value  in  all  instances.  Note  that  the 
iteration  should  not  be  started  with  the  zero  matrix  because  it  is  a  stationary  point  of 
the  sum  of  squares. 


Problems  from  Salane  [1987] 
n  m 

37.  2  16  Hanson  1 

38.  3  16  Hanson  2 
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Problems  from  McKeown  [1975a]  (also  McKeown  [1975b]) 


39a. 

n 

2 

m 

3 

McKeown  1 

0.001 

39b. 

2 

3 

McKeown  1 

0.01 

39c. 

2 

3 

McKeown  1 

0.1 

39d. 

2 

3 

McKeown  1 

1.0 

39e. 

2 

3 

McKeown  1 

10.0 

39f. 

2 

3 

McKeown  1 

100.0 

39g. 

2 

3 

McKeown  1 

1000.0 

40a.  f 

3 

4 

McKeown  2 

0.001 

40b.  f 

3 

4 

McKeown  2 

0.01 

40c.  f 

3 

4 

McKeown  2 

0.1 

40d.f 

3 

4 

McKeown  2 

1.0 

40e.f 

3 

4 

McKeown  2 

10.0 

40f.f 

3 

4 

McKeown  2 

100.0 

40g.f 

3 

4 

McKeown  2 

1000.0 

41a. 

5 

10 

McKeown  3 

0.001 

41b. 

5 

10 

McKeown  3 

0.01 

41c. 

5 

10 

McKeown  3 

0.1 

41d. 

5 

10 

McKeown  3 

1.0 

41e. 

5 

10 

McKeown  3 

10.0 

41f. 

5 

10 

McKeown  3 

100.0 

41g. 

5 

10 

McKeown  3 

1000.0 

f  In  the  data  defining  this  problem  given  in  McKeown  [1975a]  and  [1975b],  the  matrix 

/  2.95137  4.87407  -2.0506  \ 

B  =  4.87407  9.39321  -3.93181 

V -2.0506  -3.93189  2.64745  / 

is  in  error  (it  should  be  symmetric).  The  value 

/  2.95137  4.87407  -2.0506  \ 

B  =  4.87407  9.39321  -3.93189  )  , 

\ -2.0506  -3.93189  2.64745  ) 

which  is  correct  to  six  decimal  digits,  was  used  in  our  formulation  of  the  problem. 
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# 


* 


Problems  from  DeVilliers  and  Glasser  [1981]  (also  Salane  [1987]) 


n 

m 

starting  value 

42a.° 

4 

24 

DeVilliers  and  Glasser  1 

(1.0,8.0,4.0,4.412) 

42b.° 

4 

24 

DeVilliers  and  Glasser  1 

(1.0, 8.0, 8.0, 1.0) 

42c.° 

4 

24 

DeVilliers  and  Glasser  1 

(1.0,8.0,1.0,4.412) 

42d.° 

4 

24 

DeVilliers  and  Glasser  1 

(1.0, 8.0, 4.0, 1.0) 

43a. 0 

5 

16 

DeVilliers  and  Glasser  2 

(45.0,2.0,2.5,1.5,0.9) 

43b.° 

5 

16 

DeVilliers  and  Glasser  2 

(42.0,0.8,1.4,1.8,1.0) 

43c.° 

5 

16 

DeVilliers  and  Glasser  2 

(45.0,2.0,2.1,2.0,0.9) 

43d.° 

5 

16 

DeVilliers  and  Glasser  2 

(45.0,2.5,1.7,1.0,1.0) 

43e.° 

5 

16 

DeVilliers  and  Glasser  2 

(35.0,2.5,1.7,1.0,1.0) 

43f.° 

5 

16 

DeVilliers  and  Glasser  2 

(42.0,0.8,1.8,3.15,1.0) 

Problems  from  Dennis,  Gay,  and  Vu  [1985] 


44a.°  f 

n 

6 

m 

6 

Exp. 

791129 

starting  value 

(.299,  -0.273,  -.474,  .474,  -.0892,  .0882)* 

44b.°  f 

6 

6 

Exp. 

791226 

(-.3,  .3, -1.2,2.69, 1.59, -1.5) 

44c.°f 

6 

6 

Exp. 

0121a 

(-.041,  .03,  -2.565,2.565,  -.754,  ,754)| 

44d.°f 

6 

6 

Exp. 

0121b 

(-  .056,  .026,  -2.991, 2.991 ,  -  .568,  .568) 

44e.°f 

6 

6 

Exp. 

0121c 

(-.074,  .013,  -3.632,3.632,  -  .289,  .289) 

45a. 0 

8 

8 

Exp. 

791129 

(.299,  .186,  -0.273,  .0254,  -0.474,  -.0892,  .0892)* 

45b.° 

8 

8 

Exp. 

791226 

(-.3,  -  .39,  .3,  -  .344,  -1.2,2.69, 1.59,  -1.5) 

45c. 0 

8 

8 

Exp. 

0121a 

(-.041,  -.775,  .03,  -.047,  -2.565, 2.565,  -.754,  .754)* 

45d.° 

8 

8 

Exp. 

0121b 

(-.056,  -.753,  .026,  -.047,  -2.991,2.991,  -.568,  .568) 

45e.° 

8 

8 

Exp. 

0121c 

( -  .074,  -  .733,  .013,  -  .034,  -3.632, 3.632,  -  .289,  .289) 

f  Variables  X2  and  x4  (b  and  d  in  Dennis,  Gay,  and  Vu  [1985])  are  eliminated  from  the 
linear  constraints  in  order  to  get  the  6-variable  formulation  of  the  problem  (see  Dennis, 
Gay,  and  Vu  [1985]). 


I  Specification  of  some  starting  values  in  Dennis,  Gay,  and  Vu  [1985]  is  incomplete.  The 
correct  values  were  obtained  from  D.  M.  Gay  in  1986. 
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